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ABSTRACT 



This thesis contains discussion, theory and program code 
for a computational fluid dynamics (CFD) model of a wing of 
arbitrary planform. The model assumes incompressible, 
inviscid, irrotational flow. The program computes forces 
acting on the wing by modeling the flow with a set of horse 
shoe vortex elements. It models the flow over an arbitrary 
wing using two solutions. One solution is the ideal lift, 
associated with a cambered and twisted wing. The other 
solution is the additional lift associated with a flat wing. 
The program computes wing camber and twist using an elliptic 
loading distribution. The thesis includes the FORTRAN source 
code, a separable User's Manual for the VORTEX program, 
discussion of the theory applied in the model, and 
instructions for operating the program. It shows a sample 
wing planform with tabular and graphic results. 

The thesis also discusses two other CFD models based on 
circulation (F) and pressure difference (ACp). It presents 
some of the problems and solutions in grid generation. 
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I . INTRODUCTION. PURPOSE AND GOALS 



AE2035 is an introductory course in aerodynamics for 
aeronautical engineering students at Naval Postgraduate 
School, Monterey, CA . The course must provide students with 
the required knowledge for future courses in the 
Aeronautical Engineering Curriculum. The student must 
understand the techniques used in computational fluid 
dynamics (CFD) to satisfy that requirement. Powerful 
computational models are available for the working 
a e r o dy nam i c i s t , but they are seldom usable in an 
introductory course. The programming language for these 
models is usually optimized for computational efficiency 
rather than teaching, and the documentation may be poor or 
unavailable. Students need a computer program with graphic 
output which is less complicated to use as a teaching tool. 
Providing that computer program is the goal of the thesis. 

The program runs on a micro or desk top computer, in 
keeping with the emphasis on use of individual workstations 
to supplement the school's mainframe computer. The micro 
computer has better graphic output programs available, 
without the extensive programming required on the school 
mainframe. This flexibility allows the student to modify the 
output to suit his/her requirements. 
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The thesis contains five principal sections. 



The first 



section, TECHNIQUES, provides a brief synopsis of three 
different flow models. The section also provides some 
insight to certain problems that arise in CFD. 

The second section, VORTEX PROGRAM DEVELOPMENT, contains 
a detailed discussion of the horse shoe vortex model. The 
section includes the theory and techniques used in the 
VORTEX program. 

The third section, CIRCULATION MODEL, contains 
discussion of a circulation (F) model. The section includes 
the theory behind the model and some of the problems 
encountered . 

The fourth section, PRESSURE DIFFERENCE MODEL, contains 
discussion of the pressure difference (ACp) model. The 
section also includes the theory behind the model and 
problems encountered. 

The fifth and last principal section, the USERS MANUAL, 
is separable from the body of the thesis and provides 
detailed theory on the horse shoe vortex model. The section 
also includes instruction on use of the program and an 
example wing planform with the resulting data. 

The computer listing for the VORTEX program is included 
in the Appendix. 
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A. PROGRAM ASSUMPTIONS AND REQUIREMENTS 

To keep the model streamlined, the following assumptions 
were made and requirements established. 

1. Flow is steady, inviscid, i r r o t a t i o n a 1 , and 
incompressible . 

2. The wing surface has zero thickness. 

3. The planforra has any sweep, taper, and aspect ratio. 

4. The program must support low aspect ratio planforms 

(less than 1.0). 

5. The program must contain comments, for easy 

modification and change. 

6. The program results must be available in tabular and 

graphic form. 

7. The program must use accepted variable names, (aspect 

ratio, taper ratio, etc.) 

8. The program must have two options. One is to generate 

loading for a flat wing shape. The other is to 

generate a wing shape for an elliptic loading. 

9. The programming language must be FORTRAN, with 

executable and source files provided. 
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II . TECHNIQUES 



Three different models were investigated. Each of the 
models uses a governing equation for the induced velocity at 
each control point^ on the wing. That velocity is 
associated with the vortex strengths located at field points 
on the wing. The induced velocity is a function of two 
factors. One is the vortex strength and the other is the 
physical orientation of the field and control points. Flow 
must be tangent to the wing surface and is found by adding 
the induced velocity vector (w) to the remote velocity 
vector (Uoo). The end result is a direct correlation between 
the wing's shape^ and the strengths of the vortices. This 
relationship allows either the strengths of the vortices or 
the wing shape to be the independent variables. 



^Field points are the positions of the vortex elements, 
and control points are the positions where the velocity is 
evaluated. For example, a control point is where flow 
tangency is enforced. 

^Wing shape will be used throughout to describe the 
surface slope of the wing. For example, one wing shape 
might be a flat plate, another might have the tip twisted 
relative to the root chord. 
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The three models use circulation (F), 



inc r eraen t a 1 



circulation (AF), and pressure difference coefficient (ACp). 
These quantities are related as follows. 




Qt> CZ> 



A. HORSE SHOE VORTEX MODEL 

This program models the flow by a set of horse shoe 
vortices distributed over the wing. Each vortex consists of 
a bound portion perpendicular to the remote velocity, plus 
two trailing portions. Each vortex element generates an 
incremental circulation (AF) around the element. 

The program generates a matrix of influence coefficients 
from an equation for the induced velocity. The program 
finds the strength of each bound vortex and therefore the 
incremental circulation around the element by solving the 
set of simultaneous equations for wing shape. The program 
can also invert the problem and solve for the wing shape 
from a specified vortex distribution. 

The program can then find the forces and moments on the 
wing, knowing the distribution of incremental circulation. 
This solution is possible by using the Ku 1 1 a - J oukowski 
theorem which relates incremental circulation and effective 
velocity to the force generated. 
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B. CIRCULATION MODEL 



This program models the flow by a series of circulation 
elements distributed over the wing. There is also a 
circulation sheet that extends from the trailing edge to 
infinity. 

Solution of a set of governing equations provides the 
strength of each circulation element. The equations satisfy 
flow tangency on the wing surface and two boundary 
conditions. Firstly, circulation strength is zero along 
the leading edge and wing tips. Secondly, the derivative of 
circulation with respect to chord-wise direction is zero at 
the trailing edge. The derivative of circulation with 
respect to chord is vorticity. This trailing edge boundary 
condition satisfies the Kutta condition of zero vortex 
strength at the trailing edge. Zero vorticity at the 
trailing edge also means that the circulation strength 
behind the wing is a function of the span-wise coordinate 
only . 

C. PRESSURE DISTRIBUTION MODEL 

The final model of the flow uses a series of pressure 
difference elements, ACp , distributed over the wing. The 
program builds a set of N equations in N unknowns from a 
governing equation and from the requirement for flow 
tangency. The boundary conditions are: 

1. ACp is zero at the trailing edge, and along the wing 

tips. 
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2. For a flat wing, the strength of ACp along the 
leading edge is unrestrained. 

3. For a wing of elliptic lift distribution, the 
strength of ACp along the leading edge is zero. 

D. PROBLEMS INHERENT TO ALL MODELS 

1 . Field and Control Point Placement 

The placement of field and control points within 
elements is important for all models. Coincident field and 
control points exhibit singularities. Each model uses 

different methods to handle them. Direct integration 

eliminates the singularity in two models. The other 

requires a special arrangement of control and field points 
along the leading edge and wing tips. 

2 . Grid Size and Solution Speed 

Fine grids require more time for solution and larger 
computers. All models compromise between data quantity and 
speed. The circulation model is very memory intensive and 
will only run on the mainframe computer. The other models 
will run on a desktop computer, though the time required 
increases rapidly with finer grids. 

3 . Grid Element Shapes 

Finite elements in a non- rectangular planform can be 
either trapezoids, or rectangles. Trapezoids match the 
planform exactly; rectangles present a ragged leading and/or 
trailing edge. The form to use is dependent on the type of 
model. The circulation model can use a trapezoid without 
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difficulty. Coordinate transformations easily handle the 
shape and at the same time allow non-uniform cosine spacing 
of the elements. The Vortex and ACp models work best with 
rectangular grid elements. Therefore the planform is not an 
exact representation of the wing. This modification is 
reasonable since the forces and moments computed are 
comparable to other models. 

E. SCOPE OF DISCUSSION 

Of the three models investigated, only the VORTEX 
program was completely successful. Significant effort was 
expended on the CIRCULATION and ACp models trying to make 
them functional. Rather than expend a large portion of this 
thesis discussing all aspects of the two unsuccessful 
models, an overview of them will be presented and the 
relevant factors of each will be addressed. Since the 
VORTEX program was successful, and is completely functional, 
a detailed discussion of the procedure and assumptions is 
provided. In keeping with this, the FORTRAN source code for 
the VORTEX model is the only computer code included. 



8 



Ill . VORTEX PROGRAM DEVELOPMENT 



A. BASIS OF THE VORTEX MODEL 

The VORTEX program is an adaptation of a program by 
Moran [Ref,l]. The program models the flow by a series of 
horse shoe vortices distributed over the wing. In his 
text, Moran [Ref.l] develops a simple program to find the 
strengths of a series of horse shoe vortices. His program 
is rather limited in that it only works for straight wings 
without taper or sweep. Additionally, his program only 
finds the loading generated by a flat wing. The final 
limitation is his use of uniform sized grid elements, which 
fails to concentrate grid elements near the boundaries of 
the wing . 

As with many aerodynamic problems, the VORTEX program 
uses a grid or mesh of finite sized elements distributed 
over the surface of the wing. The program divides the wing 
planform (there is zero thickness modeled in this program) 
into N discrete elements. 

The program solves a set of N equations for N unknowns. 
The equations are derived from an induced velocity equation 
and the requirement for flow tangency on each grid element. 
The N unknowns are the strengths of the individual horse 
shoe vortices associated with each element. 
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Knowing the strengths of the vortex elements, the . 
program can find the forces and moments acting on the wing 
as a result of those elements. 



B. MODIFICATIONS TO THE MODEL 

The VORLAT program, developed by Moran [Ref.l], required 
changes to suit the goals of the thesis. Items requiring 
improvement were the ability to handle sweep and taper. 
Other improvements were also added to generate the grid 
using cosine spacing. A module to determine camber and 
twist, called shape, was also incorporated. 



C. TECHNIQUE FOR A FLAT WING 

This is a short description of the VORTEX program. 
Consult Moran [Ref.l] for a description of the * VORLAT 
program, the foundation of the VORTEX program. 

1 . Down-wash and the Relationship to Flow Tangencv 

A horse shoe vortex is located in the xy plane as 
shown in Figure 3-1. The vortex induces a velocity at any 
point in the xy plane. That velocity can be determined from 
equation 3.1, as follows: 



Ar 

^ (^,y) = 

ij 4n(y-y ) 



V(x-i f+(y-y 



1 + 



(x-x ) 
a 



Ar 



4nCy-3'^) 



1 + 



V(i-i r 

a 0 

ix — x) 
a 



3 . 1 
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If X 



Xq then, 



w 

u 



AP / 1 

4ii \ V— V 

a 




3 . 2 



^ij(x,y) is the velocity induced at the control point (x,y), 
AT is the incremental circulation^, and (x^,y^) and (x^,yi^) 
are the coordinates of the corners of the horse shoe vortex. 
Moran derived this equation using the Biot-Savart Law 
[Ref.l]. ^ 




A Horse Shoe Vortex in the Wing Plane 
Figure 3-1 



^Different authors use different variable names for 
Circulation and Vorticity. In the Vortex model, AT is used 
for incremental circulation. In the Circulation model, T is 
used for Circulation strength. Caution is urged when 
relating one model to another. 
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As shown in Figure 3-2, the remote velocity, with a wing 
at angle of attack a, has components VcoCOs(a) and VooSin(a). 



Vco 

Vco cos(cO FLAT VING 

Velocity Components on a Flat Wing 
F i gur e 3-2 

If the flow is tangent to the wing, the down-wash 

velocity^ at a point on the wing must equal VcoSin(a). Each 

vortex on the wing contributes to the down-wash^ at every 

point on the wing. So, an equation can be developed for each 

control point as a function of the strength of each horse 

shoe vortex. In this thesis, the mid point of any bound 

vortex is termed a field point. Any point on the wing 

\ 

where flow tangency is evaluated is termed a control point. 




9 ... 

“^Downwash is considered positive when its direction is 

downward . 

^The contribution will be positive when the induced 
velocity is downward, or negative when the induced velocity 
is upward. 
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An example of the equation for the control point (Xp,yp) is. 

V sina — ^ V w (x j )= 0 ^ ^ 

® — ij ij p'^p 

u 

Combining all the control points results in a linear set of 
N equations in N unknowns. 

2 . Placement of the Horse Shoe Vortex 

Moran [Ref.l] and others have discussed the 

placement of a lumped vortex on a two-dimensional airfoil 
with one chord-wise element. For the two-dimensional case, 
the vortex is placed at the quarter chord point and flow 
tangency is evaluated at the three-quarter chord point. The 
two-dimensional reasoning can be applied to a wing of finite 
span, assuming the wing is formed with a single choYd-wise 
element. 




Circulation Distribution over a Flat Airfoil 

Figure 3-3 
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For the flat airfoil, the circulation distribution 
is known to be approximately as shown in Figure 3-3. The 
centroid of this circulation distribution is at the quarter 
chord point. As a result, it is reasonable to lump the 
total circulation at the quarter chord point. For these 
conditions, the flow is tangent at the three quarter chord 
point. This arrangement correctly computes the moment 
coefficient (C^) for the two dimensional flat airfoil. 



The 


reason for 


other 


authors ' 


placement of 


the 


vortex at 


the quarter 


chord 


o f each 


gr id e lement 


was 


que s t i one d 


during development 


of the VORTEX program. 


The 



other authors used lifting line theory with a single chord- 
wise element for their quarter chord vortex placement. The 
VORTEX program does not use a single chord-wise element. 
The wing is divided into a number of individual elements 
which are modeled with a uniform distribution of vorticity 



over each 


element . 


The 


centre id is 


a t 


the center of 


the 


e lement . 


So , for 


the 


arrangement 


in 


F i gure 


3-4, 


the 


incremental 


c ircula t ion 


is concentrated 


a t the 


center 


o f 



each e lement . 
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Approximation of Vorticity with Individual Elements 

Figure 3-4 

With the circulation concentrated at the center of 
each element, the flow tangency point is set at one half an 
element downstream of the concentrated circulation and falls 
at the border between elements, or at the trailing edge of 
the wing. 

N 

The coefficients affected by the vortex placement 
are and X^p . In Ref. 2, X^p/c for a straight flat wing, 

aspect ratio 2, was determined to be 0.209. With the lumped 
vortex at the mid chord point, the location of X^p/c 



converged 


to a 


value 


of approximately 


0.234 


The error 


indicates 


the 


vortex 


is too far 


aft 


on 


the element. 


Relocating 


the 


lump e d 


vortex element 


to 


the 


quarter chord 


po int , and 


evaluating 


flow tangency 


a t 


the 


three - quarter 



chord point, moved the computed location of X^p/c forward, 
but still not to the true value of 0.209. 



15 



In Ref. 3, Hough looked at optimum grid and vortex 
arrangement for rapid convergence. One of the planforms 
was a straight flat wing with an aspect ratio of 2. Hough 
used a conventional vortex lattice arrangement of uniform 
size grid elements with the lumped vortex at the quarter 
chord point. For the aspect ratio 2 wing, he found that 
the computed Cj^/a was larger than actual and that the value 
converged as in Figure 3-5. It was expected that the error 
in Ci^/a could be reduced by decreasing the planform area by 
some factor. He found that convergence of Ci^/a , Vortex Drag 
Factor (K)^, and X^p was improved by insetting the tip 
vortex by some fraction (d) of a grid element. This inset 
improved convergence dramatically when d = 1/^. Though not 
addressed by Moran [Ref.l], that is the reason for insetting 
the grid by d = 1/4. 

The VORTEX program uses cosine spacing instead of 
uniform spacing so insetting the grid by d = 1/4 was not an 
available option. Instead, the span-wise grid layout is 
developed with an additional row of tip elements which are 
not used in the computation of incremental circulation or 
force. This reduces the planform area and the improvement 
in convergence can be seen in Figures 3-5 and 3-6. 



•■Vortex Drag Factor, K is; 7r*AR*Cj)j_ 



Cl' 
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Percent error in CL/o 



Additionally, the convergence of X 
Figure 3 - 7 . 



s een in 



cp is improved, 




Percent Error in CL/a 
F i gur e 3-5 
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Vortex Drag Factor (K) 
Figure 3-6 
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(Number of Chordwise vortices)” 

^cp 

Figure 3-7 

3 . Satisfaction of the Kutta Condition 

To satisfy the Kutta condition, the vortex strength 
must be zero at the trailing edge. While not explicitly 
required by the VORTEX program, the resultant vortex 
strengths, for a flat wing, approach zero at the trailing 



edge. Therefore, the results seem to imply that flow 
tangency at the trailing edge of a flat wing is a corollary 
of the Kutta condition. 

4 . Planform Development 

It is necessary to develop wing geometry from aspect 
ratio, taper ratio and sweep angle. Aspect ratio is the 
span divided by the average chord. Taper ratio is the tip 
chord divided by the root chord, and sweep angle is the 
angle between the leading edge and the y axis. The y axis 
is perpendicular to the remote velocity and parallel to the 
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earth. The aircraft design courses taught at NFS use these 
variables to define the planform. Moran [Ref.l] uses a 
fixed root chord equal to 1.0 and varies the span to achieve 
different aspect ratios. The VORTEX program adopts this 
approach, which works well. 




Planform and Variables used in VORTEX Program 

F i gur e 3-8 

The subroutine SET85 determines planform variables 

\ 

and stores them in part of the matrices WING and SECTN for 
further use. The WING and SECTN matrices also contain final 
results of the program. 

5 . Grid Development 

The program develops a grid to model the vortex 
system after establishing the outline of the planform. The 
grid elements are rectangular. The number of elements must 
be small to keep the time required for solution reasonable. 
The program concentrates grid elements in areas where the 
rate of change is most rapid, or where the values are most 
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important. Those areas are along the boundaries. An 
excellent method for concentrating the points in the areas 
desired while minimizing the total number of points is 
through cosine spacing. The program uses a cosine function 
to distribute span-wise and chord-wise grid points in a non- 
uniform fashion after a suitable coordinate transformation. 
Figure 3-9 contains a sample of the layout. Notice how the 
rectangular elements model the planform. 



Remote Velocity 

> 




Grid Model with Cosine Spacing 
Figure 3-9 

6 . Solving the Set of Equations 

The method used to solve the set of N linear 
equations is arbitrary. Moran [Ref. 1] uses a Gaussian 
algorithm for solution of an N by N matrix. The VORTEX 

program uses the same algorithm. Increasing the total 
number of grid elements beyond 100 would require 
modification of the GAUSS subroutine. 



20 



7 . Finding the Forces on the Elements 

The Kutta- Joukowski theorem states that the force 
per unit span acting on an element is. 







3 .4 



In this equation, ^eff local effective velocity at 

the center of the element. V^ff is defined in equation 3.5 
and Figure 3-10. p is the density and AT is the incremental 
circulation around the element. This circulation (AT) is 
the same AT contained in equation 3.1, and can be termed the 
incremental circulation that occurs over the element. 
Throughout this section, certain approximations and 
dimensional simplifications will be made. The firsts is the 
small angle approximation, where the sine is approximately 
equal to the angle in radians, and the cosine is 
approximately equal to 1.0. This approximation requires 
that angle of attack for the flat wing be small, generally 
less than 10 degrees. The approximation also requires that 
any wing slope on the elliptically loaded wing be small. 
This requirement is met by restricting the desired lift 
coefficient to values less than about 0.5. In the process 
of determining coefficients of lift, drag and moment, 
certain dimensional reference quantities arise. These 

quantities are density (p), remote velocity (Vqo) , planform 
area (S) , and average chord (c) . In performing dimensional 
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analysis an arbitrarily value of unity can be assigned to 
these terms . 

a. Finding ACp 

The distribution of ACp over the wing is 
desired, so the force on each element must be converted to a 
dimens ionless pressure difference coefficient. ^eff 

defined relative to the wing as; 




-> 

= V o)sa i 

(X> 



-f 



Vjsin a — IV 



\ 



J 
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Veff 



Components of Effective Velocity (^eff^ 

Figure 3-10 

It is also significant to note that AT can be written as; 



A? = Ar j 
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since the bound portion of the horse shoe vortex is aligned 



with the y axis. The result of the force cross product, 
Equation 3.4, is ; 



Ay 



= P 



w 



V sin a 




-> 

+ V cos a A r i 

OD 
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Setting p -= 1 , V’oo “ and incorporating the small angle 

approximation; 




a; - a AT / + 
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This force is related to ACp . ACp is a scalar, rather than 
a vector quantity, so the program uses the force component 
normal to the wing to compute ACp. That component is*; 



A^’ = Ay AT 
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Making the force non-dimensional gives a pressure 
difference coefficient, ACp. 



AC = 

p 



AF 






3 . 10 
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Finally, afcer setting p = 1 and Va, ■= 1, ACp for an element 
is ; 



2AF 
AC = 

^ /\rAy 



2 Ar 

Ax 
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This equation provides the ACp for each element, which can 
be plotted versus chord, for each span-wise section, 
b. Finding the lift, drag and moment 

Lift, drag and moment are also found from the 
vector force on each element. Equation 3.7 showed that the 
force on an element has components normal and tangent to the 
wing. Once more setting p = 1, and Voo = 1, the vector 

force on an element is; 






w — sin a 



jA l^Av i -I- cos a AFAjy ^ 
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These forces are oriented relative to the wing coordinate 
system, where i is parallel to the wing and k is 

perpendicular to the wing. Lift and drag are normal and 
tangent, respectively, to the remote velocity (Vco) . The 
wing is at an angle of attack (oc) to that remote velocity. 
Using this information, the program transforms the forces to 
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a coordinate system oriented to the remote velocity. For an 
individual element the lift and drag forces are; 

AL = Elemental Lift Force = (cosa)(AyAD (cosa) — (w — sin a) (AyAD (sm a) 

AD = Elemental Drag Force = (cos a) (Ay AD (sm a) + {lu — sm a) (AyAF) (cos a) 

using a small angle approximation, and discarding 

higher order terms, 3.13 

AL = Elemental Lift F orce = AyAF — ioaAyAF = AyAF 
AD = Elemental Drag Force — AyAFa + (ic — a)AyAF = to AyAF 

The total lift and total drag acting on the wing is a 
summation of the elemental lifts and drags. Using wing 
symmetry ; 



L = Total Lift = 2 V 

semi - span 3.14 

D = Total Drag - 2 ^ AD 

semi — span 

The moment about the y axis for an individual 
element, with the accepted sign convention, is; 

AM = Elemental Moment — —(normal wing force) (x position of center of element) ^ 

AM = Elemental Moment — ( — AyAF)(x ) 

cen 
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Total moment about the y axis is a summation of 
the elemental moments. Once again, using wing symmetry. 



M — Total Moment = 2 ^ AM 

semi — span 
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8 . Finding the Coefficients 

Non- dimens ional iz ing the total lift and drag gives, 



L 




C 



D 



D 
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where S = span x average chord = (b x c), p * 1 and Vqq “ 1. 

Non - d iraens iona 1 i z ing the total moment about* the y 
axis gives , 



M 

~ I I ' 3.18 

- p c S 

where c is the average chord. 

The aerodynamic center is important. X^^ is the 

chord-wise position about which the flat wing has zero 
moment, and is found from. 



QC 



— total moment 
total lift 
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This equation assumes that total moment is taken about the y 
axis, or x=0 . When the program finds the moment on the 
elliptically loaded wing, it uses the from the flat wing 
as the reference axis and finds . 

D. TECHNIQUE FOR SPECIFIED LOADING 

1 - Finding Camber and Twist from Loading 

Elliptic loading in a span-wise direction results in 
minimum induced drag, and loading is proportional to 
circulation. At sufficiently high aspect ratio, a flat wing 
with elliptic area distribution generates elliptic loading 
but weighs more than a straight wing. Wings with straight 
leading and trailing edges also cost less to manufacture 
than elliptic wings. So, the program generates a span-wise 
elliptic lift distribution by slightly twisting the wing. 
The direct relationship between wing shape (camber and 
twist) and circulation distribution makes this generation 
possible . 
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canber 




Camber and Twist on the Elliptically Loaded Wing 

Figure 3-10 

For a thin wing, elliptic load distribution in the 
chord-wise direction requires an approximately parabolic 
camber shape. The program specifies an elliptic chord-wise 
load distribution and then finds the wings' associated 
shape . 

2 . Specifying the Lift Coefficient 

Specifying the form and amplitude of the ideal lift 

distribution over t-he wing determines the form and amplitude 

\ 

of the wing shape. A load distribution of elliptic form is 
imposed in this case and the corresponding amplitude is 
fixed by the desired ideal wing lift coefficient, 

3 . Maximum Desired Lift Coefficient 

The small angle approximation breaks down when large 
lift coefficients are specified. To prevent this 

occurrence, the maximum desired lift coefficient is 
restricted to 0.5. 
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4 . Achieving the Desired Lift Coefficient 

The program scales a reference distribution of ACp , 
which is elliptic in form, to achieve the desired total lift 
coefficient for the wing. 

5 . Forces. Moments and Coefficients 

Once the ACp distribution is scaled for the 
elliptically loaded wing, the program finds forces, moments 
and coefficients in the same manner as for the flat wing. 
The only difference is in the coordinate transformation used 
to convert forces on the wing to lift and drag. The 
individual elements on the elliptically loaded wing are no 
longer at a uniform angle to the remote velocity. This 
variation in angle must be taken into account when 
performing the coordinate transformation which gives lift 
and drag forces on individual elements. 

The program finds moment coefficient about the 
aerodynamic center,' , from, 

= + 3.20 

c ^ 

where is found from the flat wing, and c is the average 

chord. pitching moment of the elliptically 

loaded wing about the y axis, and CLi is the lift 
coefficient of the elliptically loaded wing. 
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E. VALIDATION OF RESULTS 

Two sources are used to check the validity of the VORTEX 
program results. Lifting line theory provides one source of 
predicted lift and drag values for straight high aspect 
ratio wings. A continuous loading method developed by the 
National Aerospace Laboratory of the Netherlands (NLR), 
presented in Ref. 2, is used as the second source. The 
method used by NLR is a very accurate computational method 
and is considered to be exact for this comparison. 

Kuethe and Chow [Ref. 4] discuss lifting line theory for 
the case of a flat, untapered, rectangular wing. Using 
lifting line theory and the VORTEX program, lift and drag 
coefficients for identical wings were computed. A range of 
aspect rations from 0.5 through 20 were selected. Lifting 
line theory is known to be reliable at higher aspect ratios, 
but tends to lose accuracy at low aspect ratios. 

As the graphic results in Figures 3-11 and 3-12 
demonstrate, CL/q and CDi/(a)^ values obtained from lifting 
line theory and the VORTEX program converge nicely at high 
aspect ratios. At lower aspect ratio, lifting line and the 
VORTEX program show a difference in calculated value. The 
values of CL/q and CDi/(a)^ obtained by NLR [Ref. 2] agree 
closely with the VORTEX program results at aspect ratio 2. 
The NLR results validate the VORTEX program results for 
straight wings. Exact data from swept and tapered wings 
were not investigated. 
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Comparison of VORTEX Program results 
with Lifting Line Theory 




Aspect Ratio 

Comparison of Lift 
Figure 3-11 

Comparison of VORTEX Program results 
with Lifting Line Theory 




Aspect Ratio 

Comparison of Drag 
Figure 3-12 
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IV . CIRCULATION MODEL 



A. BASIS OF THE CIRCULATION MODEL 

This program models the flow with a sheet of distributed 
circulation and a continuous trailing vortex sheet behind 
the wing. Barna [Ref. 5], and also Mi Ine - Thoms on [Ref. 6, 

pages 171-177], discuss this conceptually simple model. The 
program uses the wing shape, flow tangency, boundary 
conditions, and the Biot-Savart Law to develop a set of 
equations. The program then solves the set of equations for 
the unknown circulation strengths at all grid points on the 
wing. 

B. TECHNIQUE FOR SOLUTION OF THE CIRCULATION MODEL 
1 . Induced Velocity 

The velocity induced at a point by the distributed 

\ 

circulation sheet on the wing and the trailing vortex sheet 
can be treated as the sum of the velocity induced by each. 
The velocity induced at the control point (Xq,Yq) by field 
point (x,y), located on the wing is. 



r is 
r is 
(x,y) 



IV 



wing 



+ 1 



1 



I V 4n 









dx ^ dy 
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the distance between the field and control points, and 
the strength of the circulation at the field point 
. Setting the semi-span length equal to unity, and 
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varying the chord accommodates different aspect and taper 
ratios . 




Planform and Variables used in Circulation Model 

Figure 4-1 

Circulation (F), and pressure difference coefficient 
(ACp) are related as follows; 




on 



4 . 2 



In equation 4.2, ACp is the pressure difference coefficient 
between the upper and lower surf aces , and T is the 
circulation . 
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The velocity induced at control point (Xq,Yq) by the 
segment of the trailing vortex sheet attached to the wing at 
(^t^y) f is : 



w 



trail 



+ ' 1 I r-(x -X ) 

i I 1 O 

-1 ^y-yj 




4 . 3 



r is the distance between the points, and is the strength 
of the circulation at the field point trailing 
edge . 

When the remote velocity is unity, the induced 
velocity is the same as the slope of the wing at the point 
in question, and the total slope is: 



addi t ional ly , 





wing 



wing 

V 



and 



r) 



^ trail 



trail 

V 
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4 . 5 



4 , 6 



The program uses these equations to satisfy flow 
tangency on the wing with a known shape and solves for the 
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unknown circulation strengths. In matrix notation the 

equations have the following form. 





\ 


B 


f 


bz 1 


A 


— + 


1 — 1 “ 






l dx 1 




1 by 1 


bx 1 



2 . Boundary Conditions 

The requirement for flow tangency at each element on 
the wing is not enough to find the circulation which 
satisfies the Kutta condition; the program must also enforce 
certain boundary conditions. 

The flat wing, which produces the additional lift, 
has the following boundary conditions: 

1. r « 0 along the leading edge and the wing tips. 

2. dr/dx * 0 along the trailing edge. 

The elliptically loaded wing, which produces the 
ideal lift, has the following boundary conditions: 

1. r »= 0 along the leading edge, and wing tips. 

2. dr/dx = 0 along both the leading and trailing edges. 

Using a finite difference scheme, it is possible to 
approximate the partial derivatives of T (dT/dx, & dT/dy) 
from the values of F at neighboring elements and satisfy 
these boundary conditions. The partial derivative of F with 
respect to x, namely (dF/dx), approaches infinity at the 
leading edge of the flat wing. This singularity makes 
finite differencing near the leading edge difficult. 

Cosine spacing is used to minimize that problem and also to 
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spacing is 


used to minimize 


that 


problem 


and 


also to 


concentrate 


elements near the 


e dge s 


0 f the 


wing . 


The X 



coordinate in the chord-wise direction can be expressed in 
an alternate form by (^ , where: 



(f) = cos 



-1 



1 - 



2(x — Aoy) 
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The additional variables are defined in equation 4.24 and 
Figure 4-3. The y coordinate in the span-wise direction 
becomes 6, where: 



0 = C05 ^(y) 
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After this transformation,; 



— 

a4> 



dx 



c 

— sm(d)) 
2 



4. 10 



From the transformation, it is possible to show that dT/d<f> 
has a finite value at the leading edge even when dT/dx is 
infinite . 
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Circulation vs Phi 
Figure 4-2 

Using Figure 4-2 as a guide, the program estimates 
the value of dT/d<t> from the values of F on neighboring 
elements. The following shows the procedure for finding 
dV/d<j> at the leading edge element as a function of F on 
neighboring elements. Other points are derived in a similar 
fashion . 



dV , 

— - A, + 2A^4> + SA^ct)- 



For the first element, 



For the second element, 



For the third element. 



A4> 

3A4> 



5A4> 
<^> 3 = — 



4.11 



4 .12 
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4.13a 




wh i ch be come s 



0.5 


0.25 


0.125' 




A 






1.5 


2.25 


3.375 


- 




- 


K. 

e* 

r 


2.5 


6.25 


15.625 








I’ 



and 



A A,p = 3.75 r -.833 T + .15 P 

4>2 4>3 



A A4>“=-4.0r + 2 . or - 4 F 

4 >, -l>2 4.3 



A A4>^= ir _ ,6667 r + .2T 
'*’1 4-., 4>, 



4 . 13b 



4.13c 



4 . 14 



4 . 15a 



4 . 15b 



4 ,15c 
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Ad) / Ad) 

+ 2A,— + 3A j — 
2 2 2 



2 



when 4.16 is combined with 4.15, it becomes; 



\ cf4) / j A<J) 

The uncertainty in this approximation is of order . 

interior elements, 



.5 r + .6667 I' - .1 r 
4>, 4-2 4-3 




1 

A4) 



-.5 r. 



+ .5 r. 



h-i 



h + 1 1 



and for the trailing edge element. 




1 

A4) 



-.1522 r +.9565 r 



,8043 \\ 
4*, 



Similar relations hold in the span-wise direction, 
matrix form: 



CC 



— 1 

dcf) 1 



DD 



r 




4.16 



4.17 



For 



4 .18 



4 .19 



I n 



4.20 



4.21 
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3 . Sweep and Taper 



The program must model swept and tapered wings. It 
makes a series of coordinate transformations to accomplish 
this. The final result is an equation for the slope due to 
the circulation on the wing. 



dz 

ax 



wing 



1 

4n 



^ \ \ j y y'o ^csfn4> / 3r \ 



sinOl — 1 — 

3i}) 



2 V 30 / 
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dQ 



and for the contribution of the trailing vortex filament, 



dz 


- 1 


[" 1 


r_ - (r - X ) 

1 'p U 


dx 
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dQ 
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where 



x“ = (X - pj) 

X = (x - ) 

o o o 

p = (A - C lO 

m 

A = (a n(^) , tangent of leading edge sweep 
c = root chord 

m 

X = 1 — taper ratio 




c 



Xl = 

C = C (1 - 10 v) 
m 

o = + 1 on right wing semi — span 
o = + 1 on left wing semi — span 
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4 . Field and Control Point Placement 

Field and control points are located at the center 
of the grid elements. When each grid element contains both 
a field and control point, the circulation that results from 
solving the equations oscillates wildly and bears no 
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resemblance to the expected solution. After considerable 
probing, a modified arrangement was discovered which 
generates a reasonable function. Distribution of a "curb” 
of elements containing only field points along the tips and 
leading edge gives a satisfactory solution. The remainder 
of the grid elements contain field and control points. The 
final configuration is. 



Field and Control Points 



Field Points Only 



Modified Grid Geometry 
Figure 4-3 

This arrangement gives up some degrees of freedom, but that 
can be compensated by adding an extra row of grid elements 
along the chord and an extra column of grid elements along 
the span. 

5 . Influence Coefficient Matrices 

Equations 4.22 and 4.23 can be divided into multiple 
integrals. The partial derivatives {3T/d4> & dV/dd) are 
evaluated at the centers of the elements, and the integrals 
in 4.22 and 4.23 can be represented as a summation over the 
elements. Using these assumptions and equations 4,19 
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through 4.24 it is possible to write the set of equations in 
matrix form. 

The program builds the final matrix of influence 
coefficients from a number of subsidiary matrices. When 
equation 4.20 and 4.21 are combined with equation 4.22 and 
the cosine coordinate transformation we get. 



M 
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Once [M] is generated, it is possible to obtain the 
circulation solution vector (F). 

6 . Solving the Set of Equations 

The program uses the LEQIF subroutine from the IMSL 
FORTRAN library as a linear equation solver. 

7 . Finding the Lift. Drag, and Moment 

From the distribution of circulation strength it is 
possible to find the lift and drag. Barna [Ref. 5] and 
Mi Ine - Thoms on [Ref. 6] derive the equations for lift and 
drag, as a function of circulation and down-wash. The 
program uses these equations to get lift and drag from the 
circulation on the trailing edge elements. The lift on a 
span-wise section which lies between y and (y+Ay) is, 



AL = pVFAy 



4.26 
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and the induced drag on a section is. 



AD^ = pVTu;Ajy 4.27 

r is the circulation along the trailing edge and w is the 
down-wash at the trailing edge. The procedure for finding 
moment and center of pressure is not as simple. The program 
would need to resolve forces on the individual elements to 
find moment coefficient and center of pressure. The 
CIRCULATION program was abandoned before Fortran code to 
generate moment coefficient and center of pressure could be 
written. 

C. PROBLEMS ENCOUNTERED 

1 . Coincident vs Non - c o i nc i den t Points 

The manner used to handle s ingular i t i e s is one of 
the biggest problems faced when using finite elements to 
model continuous functions. The circulation model behaves 
very well when the field and control elements are not 
c o i nc i den t . 

Field elements induce velocities which increase very 
rapidly as the distance to the control point decreases. 
When the field and control points are in the same element it 
appears the induced velocity is indeterminate. However, 
Mi Ine - Thoms on [Ref. 6] shows that an element induces zero 
velocity on itself. This result of the singularity is 
important and requires special attention. 
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2 . Field and Control Points in the Wing Interior 

The direction of velocity induced at a control point 
by its neighboring field points is significant. Using a 
simple 3X3 element grid near the center of the wing as an 
example, Figure 4-4 shows the direction of the induced 
velocity from 8 field points. All eight field points have 
positive circulation and surround the center control point. 



REMOTE VELOCITY 
> 



down 


up 


up 


down 


B 




up 


down 


up 


up 




con*fcrol poln“fc 



Interior Control Point 
Figure 4-4 

The three elements upstream from the central control point 
induce a down-wash on the control point. The other five 
elements induce an up-wash on the control point. The net 
velocity is a sum of all eight elements. Recall that the 
velocity induced is dependant on 3F/3x, not on F itself^. 
This net velocity will be down-wash if 5F/3x in the three 
upstream elements is sufficiently greater than in the other 
five e 1 emen t s . 
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The derivative of circulation is proportional 



to ACp . 




3 . Field and Control Points at the Leading Edge 

A problem becomes apparent when checking the induced 
velocity at an element on the leading edge. For a 2X3 grid 
on the leading edge, Figure 4-5 shows the direction of 
velocity induced by each of the neighboring field points on 
the central control point. 



REMOTE VELOCITY 
> 



Leading Edge Control Point 
Figure 4-5 

For this situation, all the induced velocity is up-wash, 
assuming positive circulation. However, the wing requires a 
net down-wash for flow tangency. Using this arrangement of 
points, the solution vector oscillates wildly between 
positive and negative values. This oscillation is a result 
of the system trying to generate down-wash at the leading 
edge e laments . 

The program uses a simple solution to this problem. 
By ignoring the requirement for flow tangency on the row of 
elements along the tips and leading edge, the program gives 
a net down-wash at all other elements. This simplification 
results in a well behaved, positive circulation over the 
entire wing . 
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It appears the set of equations has more unknowns 
than equations, but boundary conditions provide the 
additional constraints. 

^ - Model Complexity 

Developing a computer program usable as a teaching 
tool for basic aerodynamics is the principal goal of the 
thesis. The CIRCULATION program does not fully achieve that 
goal. The concept is relatively understandable, but is 
difficult to execute. The complex method needed to build 
the influence coefficient matrix reduces its usefulness as a 
teaching tool. The student needs a simpler tool. The 
VORTEX program discussed in Section III is that tool. 



V . PRESSURE DIFFERENCE MODEL 



A. BASIS OF THE ACp MODEL 

Like the vortex and circulation models, the pressure 
difference model (ACp) also uses a governing integral 
equation as its foundation. The integral equation relates 
wing slope at a specific point on the wing to the 
distribution of pressure difference over the whole wing. 
Dividing the wing into N grid elements, with known slopes, 
the program solves for N pressure differences (ACp). 

The requirements of the Kutta condition are not 
explicitly satisfied by the equation. Instead, additional 
constraints are necessary at the wing tips and trailing 
edge. The simplest form of satisfaction is to require that 
ACp be zero at those grid elements. That constraint is more 
severe than necessary. ACp must be zero at the edge of the 
element but the value at the center of the edge elements can 
be finite and is approximated by fitting a polynomial 
through the edge and neighboring control points. This 
approach is similar to that used in the CIRCULATION model. 
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n 



B. TECHNIQUE FOR SOLUTION OF THE ACp MODEL 
1 . The Singularity 

The integral equation governing the wing 
slope/pressure difference relationship is, 



Z\x,y) = 






5 . 1 



where the slope, 



and , 



Z\x,y) = 



£ 

dx 



ZJxj) 



5 . 2 



k 



1 1 
8n 



(x-x,) 

1 + 

V 4 - iy-y^^ 
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When the field and control points are coincident (r=0) or 
share the same y value (y^yl)t there is a strong 
singularity. It is possible to evaluate this integral and 
eliminate the singularity. 
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2 . Evaluation of the Integral 

Using the mirror image element on the opposite semi- 
span, it is possible to reduce the integral to, 



Z“(x,y)= j AC^(Xj, 3 'j)d!xj dyj 



5.4 



where , 



and 



, <yo , , o 

k = k.+ k 



5 . 5 



> 






(x-x, 



1 + 



V(x-Xjf + 
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The singularity was eliminated by evaluating the integral. 
Using Figure 5-1 as a guide, the result is. 
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Variables used for a Grid Element 
F i gur e 5-1 
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and k 



^<^ 9 +^ 10 +^11 + ^ 12 ^ 



Q, = (x-x,) + 5 
Q^= (x- i|) - 5 
Q3 = <9' - 9',) + c 



Q 4 = Cy - 9']) - c 



Q,= VqJ + Q3^ 



Qg= \/q2 + q 



Q, = Vq| + q2 



Q« = Vq^ + q^ 



4l8 
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Q. 
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^3<?4 



- 1 



Q. 
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+ Q7 



Qo 



, the mirror image is; 



*‘ = i; '«9 + «I0 + «„ + «, 2> 
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where , 



^3 = (y + 3',) + c 

^4 = (y + y,) - c 



and all other values are the same as for k. The end result 
is , 



r = k-^k^ 



and 



S/2 



In matrix form, the equations can be written. 



r 




AC 


, 


Z° 






p 







A computer program will rapidly calculate all of these 
values . 

3 . The Kutta Condition 

Solution of the equations does not guarantee a 
result which satisfies the Kutta condition^. To enforce the 
Kutta condition, the program sets ACp at elements along the 



^ It should be noted that the requirement for flow 
tangency at the elements near the trailing edge was 
sufficient to satisfy the Kutta condition in the VORTEX 
program. However that result was discovered after the 
PRESSURE program was abandoned and no effort was made to 
extend that reasoning to the PRESSURE program. 
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wing tip and trailing edge (the nth elements) equal to a 
fraction of the value of the [n-l]th elements. The program 
enforces the condition by fitting a polynomial through the 
two points. See Figure 5-2. 



REMOTE VELOCITY 
> 



N and N-1^^ Elements 
F i gur e 5-2 

4 . Grid Spacing 

The program uses a non-uniform grid spacing, based 
on a cosine function. The elements near the wing tip are 
small. This spacing has the advantages noted earlier. 

C. PROBLEMS ENCOUNTERED 

In matrix form, the model has the form. 

r 

/c®® is the influence coefficient matrix, and z® is the wing 
slope. The model specifies some values of ACp on the left 
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side of the equality and some values of wing slope (z®) on 
the right side of the equality. This complexity prevents 
use of the normal matrix solvers. As a result, the solution 
requires a complex matrix manipulation which is not readily 
understandable, or desirable for a teaching tool. 

D. CONCLUSIONS 

The evaluation of the integral over the element has 
definite advantages for coincident field and control points. 
The evaluation eliminates a strong singularity and should 
give a result using finite elements that is very close to 
the actual property. 

The Kutta condition and method used to enforce 
satisfaction is not optimal. Further studies are needed. 



n 
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B. INTRODUCTION 



LCDR Chris L. HOLM wrote this program, as a tool for use 
in AE2035. The program is partial satisfaction of the 
requirements for the degree Aeronautical Engineer from the 
Naval Postgraduate School. The objective of the program is 
to provide a simple computer simulation of flow over a thin 
wing with low aspect ratio. High aspect ratio wings are 
better served by a lifting line model. Many advanced 
programs using elaborate methods to model viscous effects, 
compressibility, boundary layer growth, and shocks, are 
available. None of them gives a good introduction to the 
field of computational aerodynamics. This program should 
fill that need at the Naval Postgraduate School. 

If you wish to skip the theory behind the program, go to 
section I for instruction on running the program. 

There are a number of ways to model the source of 
forces acting on a body surrounded by fluid in motion. 
These include potential functions, vortex distributions, 
circulation distributions, and pressure differential 
distributions. Each model is related to the other, and each 
has advantages and disadvantages. This program uses a set 
of horse shoe vortices to approximate the flow over a wing. 
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C. SYSTEM REQUIREMENTS 



System type 
Memory 

Math Co-processor (8087) 
Graphics card 
Graphics program 



IBM PC or compatible 
256K, minimum 
Recommended, not 
required . 

Recommended, not 
required . 

X, Y graphing program 
like GRAPHER, is 
recommended, not 
required . 



D. CONCEPT DESCRIPTION 

1 . The Horse Shoe Vortex 

This program, which models the flow by a series of 
horse shoe vortices, distributed over the wing, is an 
adaptation of the VORLAT program by Moran [Ref.l]. He 
develops a program to find the strengths of a series of 
horse shoe vortices associated with a flat rectangular wing. 
The VORLAT program has its foundation in two-dimensional 
airfoil theory, and Moran [Ref.l] adapted the theory to 
wings of finite aspect ratio. 

The Users Manual will present a short description of 
the VORTEX program. For coverage of the VORLAT program, the 
foundation of the VORTEX program, consult Moran [Ref.l]. 
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a. Down-wash and the Relationship to Flow Tangency 
The velocity induced at a point in the xy plane 
by a horse shoe vortex, also located in the xy plane, can be 
determined from; 

y/(x-x/+(y-3// j AF ( V(x-x/ 

(x-x ) J 4n(y— yj I (x — x ) J 

a ^ 0 a 



AF 

w.ixj) = 

ij 4n(v-v 



if ('x - x^) , 



w. ixj] 

ij 



AT / 1 

4n 






7 . 2 



w is the velocity induced at (x,y), called a control point. 
AT is the strength of the horse shoe vortex, and and 

(^a*7b) corners of the horse shoe vortex. The 

center of the horse shoe vortex is a field point. 

The velocity at a point on a flat wing, which 
is at an angle of attack a, has components V^cos(a) and 
V^s in (a) as indicated in Figure 7-1. 




Velocity Components on a Flat Wing 
Figure 7-1 
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In order to ensure the flow is tangent to the 
wing, the induced velocity or down-wash^ at that point on 
the wing must also be V^sin(a). Each vortex on the wing 
contributes to the down-wash"^ at every control point on the 
wing. So, an equation for each control point, as a function 
of all the field points can be written. An example of the 
equation for the control point (Xp,yp) is. 

Vjina ^ ^ w ix ,y ) = 0 ^ . 

® — u u P '^p 1 . 3 

y 

When the control points are combined, the result is a set of 
N equations in N unknowns. 

b. Placement of the Horse Shoe Vortex 

The bound portion of the vortex (that portion 
perpendicular to the onset flow) is placed at the 1/4 chord 
point of each grid element on the wing. Flow tangency is 
evaluated at the 3/4 chord point of each grid element. 



^Down-wash is considered positive when its direction is 
downward . 

^The contribution will be positive when the induced 
velocity is downward, or negative when the induced velocity 
is upward. 
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c. Satisfaction of the Kutta Condition 

To satisfy the Kutta condition, the vortex 

strength must go to zero at the trailing edge. For wings 
and airfoils of zero thickness, flow tangency enforced at 
the trailing edge appears to be a corollary of the Kutta 
condi t ion . 

2 . Developing the Wing 

a. Planform Geometry 

It is necessary to develop wing geometry from 

aspect ratio, taper ratio and sweep angle. Aspect ratio is 
the span divided by the average chord. Taper ratio is the 
tip chord divided by the root chord, and sweep angle is the 
angle between the leading edge and the y axis. Moran 

[Ref.l] uses a fixed root chord and varies the span as 

necessary to achieve the required aspect ratio and taper 
ratio. The VORTEX program uses this approach, which works 
well. 
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Planform and Variables used in VORTEX Program 

Figu.re 7-2 



Aspect 


ratio , AR , 


is (b^/S) 


. Span 


(b) 


i s 


the 


distance from 


wing tip to 


wing tip, 


and area 


(S) 


i s 


the 


total planform 


area of the 


wing . 










Taper 


ratio , A , is 


the ratio 


of tip chord 


to 


root 



chord . 



^ _ tip chord 7 4 

root chord c 

r 

Sweep angle delta, A, is the angle that the leading 
edge makes with a line drawn perpendicular to the root 
chord . 

b. Grid Geometry 

Using aspect ratio, taper ratio, and sweep 
angle, the subroutine SET85 develops a grid. The program 
stores the coordinates of all necessary points in the 
matrices WING and SECTN. Results of the program are also 
stored in WING and SECTN which are written to data files 
when the program finishes. 
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In a manner analogous to that used by Hough [Ref. 3], 
the tip vortices are inset to improve accuracy of the 
results. The inset distance is one element wide. 

The model keeps matrix size and the number of grid 
elements small to keep the time necessary for solution 
within reason. The model concentrates grid elements in 
areas where the rate of change is rapid, or where the values 
are most important. Those areas are along the wing 
boundaries. An excellent method for concentrating the 
points, while keeping the total number of points at a 
minimum, is through cosine spacing. Figure 7-3 contains a 
sample of the layout. Notice that the rectangular elements 
do not model the planform exactly, but in the limit the 
difference will be negligible. 



Remote Velocity 

> 




Grid Model with Cosine Spacing 
F igure 7-3 
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3 . 



Matrices Used in VORTEX Program 



The program uses two identical matrices of 
influence coefficients, [A] and [AA]. [A] is used to solve 
for the AT of the flat wing, and [AA] is used to solve for 
the wing shape of the elliptically loaded wing.^ 



A 





(FlatV/ing) 



1 . 5 



AA 




dz 1 

— } (Elliptically Loaded V^ing) 
dx 1 



7 . 6 



The SECTN matrix contains values for the wing 
sections, (chord, and width) as well as some of the final 
results. Each row of the SECTN matrix corresponds to a 
wing section in the semi-span. The Output Variables 

section (H.l.a) describes the columns. 

The WING matrix contains wing coordinates required 
by the program as well as some of the final results. Each 
row of the matrix corresponds to an element in the wing 
semi-span. The Output Variables section (H.l.b) describes 
the c o lumns . 



^The method used to develop wing shape will be 
discussed in section VII. E. 
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4. 



Solution of the Set of Equations 



The subroutine GAUSS solves the N equations in N 

unknowns. The subroutine destroys [A] in the process but 

returns the solution vector (AT) in place of the right hand 

sides. 

5 . Finding the Forces on the Elements 

The Kut t a - J oukowski theorem states that the force 

per unit span acting on an element is. 




eff 



■ y 
X Ar 



7 . 7 



In this equation, local effective velocity at 
the center of the element. ^eff defined in equation 7.8 
and Figure 7-4. p is the density and AT is the incremental 
circulation around the element. This circulation (AT) is 
the same AT contained in equation 7.1. AT can be termed the 
incremental circulation that occurs over the element. 

Throughout this section, certain approximations and 
dimensional s implif icat ions will be made. The first is the 
small angle approximation, where the sine is approximately 
equal to the angle in radians, and the cosine is 
approximately equal to 1.0. This approximation requires 
that angle of attack for the flat wing be small, generally 
less than 10 degrees. The approximation also requires that 
any wing slope on the elliptically loaded wing be small. 
This requirement is met by restricting the desired lift 
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coefficient to values less than about 0.5. In the process 
of determining coefficients of lift, drag and moment, 
certain dimensional reference quantities arise. These 

quantities are density (p), remote velocity (Voo) , planform 
area (S), and average chord (c). In performing dimensional 
analysis an arbitrarily unit value can be assigned to these 
terms . 

a. Finding ACp 

The distribution of ACp over the wing is 
desired, so the force on each element must be converted to a 
dimensionless pressure difference coefficient. ^eff 

de fined as . 



V = V cos a / + V sin a - a; I / 7.8 

eff ^ \ ^ 




Veff 



Components of Effective Velocity 
Figure 7-4 
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It is also significant to note that AT can be written as, 



^ 

AP = M j 



1 . 9 



since the bound portion of the horse shoe vortex is aligned 
with the y axis. The result of the force cross product, 
Equation 7.7, is. 



Ay 



w 



V sin a 

CD 




+ V cos a A r k 

CD 



1 . 10 



Setting p ^ 1, Voo = 1, and incorporating the small angle 

approximation. 



/ \ 4 A 7.11 
( — a jAP i -I- A P 

This force is related to ACp . ACp is a scalar, rather than 
a vector quantity, so the program uses the force component 
normal to the wing to compute ACp. That component is. 




AF = AyAP 



7 .12 



Making the force non-dimensional gives a pressure 
difference coefficient, ACp. 



AC = 

p 



AF 



-pV^^AxAy 



7 .13 
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Finally, after setting p 
is . 



1 and Voo =* 1, ACp for an element 



AC = — = — 7.14 

P AxAy Ax 

This provides the ACp for each element, which can be plotted 
versus chord, for each span-wise section. 

b. Finding the Lift, Drag and Moment 

Lift, drag and moment are also found from the 
vector force on each element. Equation 7.10 showed that the 
force on an element has components normal and tangent to the 
wing. Once more setting p = 1, and « 1, the vector 

force on an element is. 



^F = 



w— smajAPAyi +cx)saArAy/f 



7 .15 



These forces are oriented relative to the wing coordinate 
system, where i is parallel to the wing and k is 
perpendicular to the wing. Lift and drag are normal and 
tangent, respectively, to the remote velocity (Vqo) . The 
wing is at an angle of attack (oc) to that remote velocity. 
Using this information, the program transforms the forces to 
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wing tip and trailing edge (the nth elements) equal to a 
fraction of the value of the [n-l]th elements. The program 
enforces the condition by fitting a polynomial through the 
two points. See Figure 5-2. 



REMOTE VELOCITY 
> 



N and N-1^^ Elements 
F i gur e 5-2 

4 . Grid Spacing 

The program uses a non-uniform grid spacing, based 
on a cosine function. The elements near the wing tip are 
small. This spacing has the advantages noted earlier. 

C. PROBLEMS ENCOUNTERED 

In matrix form, the model has the form. 




/c®® is the influence coefficient matrix, and z® is the wing 
slope. The model specifies some values of ACp on the left 
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side of the equality and some values of wing slope (z®) on 
the right side of the equality. This complexity prevents 
use of the normal matrix solvers. As a result, the solution 
requires a complex matrix manipulation which is not readily 
understandable, or desirable for a teaching tool. 

D. CONCLUSIONS 

The evaluation of the integral over the element has 
definite advantages for coincident field and control points. 
The evaluation eliminates a strong singularity and should 
give a result using finite elements that is very close to 
the actual property. 

The Kutta condition and method used to enforce 
satisfaction is not optimal. Further studies are needed. 



54 
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B. INTRODUCTION 



LCDR Chris L. HOLM wrote this program, as a tool for use 
in AE2035. The program is partial satisfaction of the 
requirements for the degree Aeronautical Engineer from the 
Naval Postgraduate School. The objective of the program is 
to provide a simple computer simulation of flow over a thin 
wing with low aspect ratio. High aspect ratio wings are 
better served by a lifting line model. Many advanced 
programs using elaborate methods to model viscous effects, 
compressibility, boundary layer growth, and shocks, are 
available. None of them gives a good introduction to the 
field of computational aerodynamics. This program should 
fill that need at the Naval Postgraduate School. 

If you wish to skip the theory behind the program, go to 
section I for instruction on running the program. 

There are a number of ways to model the source of 
forces acting on a body surrounded by fluid in motion. 
These include potential functions, vortex distributions, 
circulation distributions, and pressure differential 
distributions. Each model is related to the other, and each 
has advantages and disadvantages. This program uses a set 
of horse shoe vortices to approximate the flow over a wing. 
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SYSTEM REQUIREMENTS 



System type 
Memory 

Math Co-processor (8087) 
Graphics card 
Graphics program 



IBM PC or compatible 
256K, minimum 
Recommended, not 
required . 

Recommended, not 
required . 

X, Y graphing program 
like GRAPHER, is 
recommended, not 
required . 



D. CONCEPT DESCRIPTION 

1 . The Horse Shoe Vortex 



This 


program , 


which models 


the flow by a series 


o f 


horse shoe 


vortices , 


dis tr ibuted 


over the 


wing, is 


an 


adap t a t ion 


of the VORLAT program 


by Moran 


[Ref.l] . 


He 


develops a 


program to 


find the strengths of 


a series 


o f 



horse shoe vortices associated with a flat rectangular wing. 
The VORLAT program has its foundation in two-dimensional 
airfoil theory, and Moran [Ref.l] adapted the theory to 
wings of finite aspect ratio. 

The Users Manual will present a short description of 
the VORTEX program. For coverage of the VORLAT program, the 
foundation of the VORTEX program, consult Moran [Ref.l]. 
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a. Down-wash and the Relationship to Flow Tangency 
The velocity induced at a point in the xy plane 
by a horse shoe vortex, also located in the xy plane, can be 
determined from: 

V(x-xf+(y-yf . , V(x-xf +(y-yf . ^ 

1 + - 1 + 

ix-x ) J 4n(y-y.) I {x-x ) J 

a 0 a 



Ar 

w (x,y) = 

y •' 4n(v-v 



if (x - x^) , 



wAx.y) = 



AT 

4n 



1 






7 . 2 



w is the velocity induced at (x,y), called a control point. 
AT is the strength of the horse shoe vortex, and and 

(^a>7b) corners of the horse shoe vortex. The 

center of the horse shoe vortex is a field point. 

The velocity at a point on a flat wing, which 
is at an angle of attack a, has components V^cos(a) and 
V^s in (a) as indicated in Figure 7-1. 




Velocity Components on a Flat Wing 
Figure 7-1 
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In order to ensure the flow is tangent to the 
wing, the induced velocity or down-wash^ at that point on 
the wing must also be VooSin(a), Each vortex on the wing 
contributes to the down-wash^ at every control point on the 
wing. So, an equation for each control point, as a function 
of all the field points can be written. An example of the 
equation for the control point (Xp,yp) is. 



Vjina - ^ AT IV ix ,y ) = 0 ^ ^ 

“ — u u P '^P 1 . 3 

When the control points are combined, the result is a set of 
N equations in N unknowns. 

b. Placement of the Horse Shoe Vortex 

The bound portion of the vortex (that portion 
perpendicular to the onset flow) is placed at the 1/4 chord 
point of each grid element on the wing. Flow tangency is 
evaluated at the 3/4 chord point of each grid element. 



^Down-wash is considered positive when its direction is 
downward . 

^The contribution will be positive when the induced 
velocity is downward, or negative when the induced velocity 
is upward. 
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c. Satisfaction of the Kutta Condition 

To satisfy the Kutta condition, the vortex 
strength must go to zero at the trailing edge. For wings 
and airfoils of zero thickness, flow tangency enforced at 
the trailing edge appears to be a corollary of the Kutta 
c ondi t ion . 

2 . Developing the Wing 

a. Planform Geometry 

It is necessary to develop wing geometry from 
aspect ratio, taper ratio and sweep angle. Aspect ratio is 
the span divided by the average chord. Taper ratio is the 
tip chord divided by the root chord, and sweep angle is the 
angle between the leading edge and the y axis. Moran 

[Ref.l] uses a fixed root chord and varies the span as 
necessary to achieve the required aspect ratio and taper 
ratio. The VORTEX program uses this approach, which works 
we 1 1 . 
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Renote Velocli:y 

> 



Planform and Variables used in VORTEX Program 

Figure 7-2 



Aspect 


ratio , AR , 


is (b^/S) 
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and area 
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area of the 
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Taper 


ratio , A , is 


the ratio 


of tip chord 


to 


root 



chord. 

^ _ lip chord ^ 7 4 

rod chord c 

r 

Sweep angle delta, A, is the angle that the leading 
edge makes with a line drawn perpendicular to the root 
chord . 

b. Grid Geometry 

Using aspect ratio, taper ratio, and sweep 
angle, the subroutine SET85 develops a grid. The program 
stores the coordinates of all necessary points in the 
matrices WING and SECTN. Results of the program are also 
stored in WING and SECTN which are written to data files 
when the program finishes. 



ct 
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In a manner analogous to that used by Hough [Ref. 3], 
the tip vortices are inset to improve accuracy of the 
results. The inset distance is one element wide. 

The model keeps matrix size and the number of grid 
elements small to keep the time necessary for solution 
within reason. The model concentrates grid elements in 
areas where the rate of change is rapid, or where the values 
are most important. Those areas are along the wing 
boundaries. An excellent method for concentrating the 
points, while keeping the total number of points at a 
minimum, is through cosine spacing. Figure 7-3 contains a 
sample of the layout. Notice that the rectangular elements 
do not model the planform exactly, but in the limit the 
difference will be negligible. 



Remote Velocity 

> 




Grid Model with Cosine Spacing 
Figure 7-3 
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3 . Matrices Used in VORTEX Program 

The program uses two identical matrices of 
influence coefficients, [A] and [ AA ] . [A] is used to solve 
for the AT of the flat wing, and [AA] is used to solve for 
the wing shape of the elliptically loaded wing.^ 



A 





{Flat Wing) 



1 . 5 



AA 




— ) {Elliptically Loaded Wing) 
dx \ 



1 . 6 



The SECTN matrix contains values for the wing 
sections, (chord, and width) as well as some of the final 
results. Each row of the SECTN matrix corresponds to a 
wing section in the semi-span. The Output Variables 

section (H.l.a) describes the columns. 

The WING matrix contains wing coordinates required 
by the program as well as some of the final results. Each 

row of the matrix corresponds to an element in the wing 
semi-span. The Output Variables section (H.l.b) describes 
the c o lumns . 



^The method used to develop wing shape will be 
discussed in section VII. E. 
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4. 



Solution of the Set of Equations 



The subroutine GAUSS solves the N equations in N 

unknowns. The subroutine destroys [A] in the process but 

returns the solution vector (AT) in place of the right hand 

sides . 

5 . Finding the Forces on the Elements 

The Ku t ta - J oukowski theorem states that the force 

per unit span acting on an element is. 






Ay 




X Ar 



7 . 7 



In this equation, ^ eff local effective velocity at 
the center of the element. ^eff defined in equation 7.8 
and Figure 7-4. p is the density and AT is the incremental 
circulation around the element. This circulation (AT) is 
the same AT contained in equation 7.1. AT can be termed the 
incremental circulation that occurs over the element. 

Throughout this section, certain approximations and 
dimensional simplifications will be made. The first is the 
small angle approximation, where the sine is approximately 
equal to the angle in radians, and the cosine is 
approximately equal to 1.0. This approximation requires 
that angle of attack for the flat wing be small, generally 
less than 10 degrees. The approximation also requires that 
any wing slope on the elliptically loaded wing be small. 
This requirement is met by restricting the desired lift 
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coefficient to values less than about 0.5. In the process 
of determining coefficients of lift, drag and moment, 
certain dimensional reference quantities arise. These 

quantities are density (p), remote velocity (Vqq), planform 
area (S), and average chord (c). In performing dimensional 
analysis an arbitrarily unit value can be assigned to these 
terms . 

a. Finding ACp 

The distribution of ACp over the wing is 
desired, so the force on each element must be converted to a 
dimensionless pressure difference coefficient. ^eff 

defined as . 



V „ = V 005Q / + V a — u; / 7.8 

eff ^ \ CO / ^ 



Veff 

Components of Effective Velocity (^eff^ 

Figure 7-4 
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It is also significant to note that AT can be written as, 



AI^ = AP 1 



7 . 9 



since the bound portion of the horse shoe vortex is aligned 
with the y axis. The result of the force cross product, 
Equat ion 7.7, is . 



Aj 



= P 




V sin a 
00 




-f V cos a A r k 

CD 



1 . 10 



Setting p = 1, V’oo = a.nd incorporating the small angle 

approximation. 



/ \ 4 7.11 
( — Q j A r / + A r 

This force is related to ACp . ACp is a scalar, rather than 
a vector quantity, so the program uses the force component 
normal to the wing to compute ACp. That component is. 




AF = AyAF 



7 .12 



Making the force non-dimensional gives a pressure 
difference coefficient, ACp. 



AC = 

p 



AF 



-pK^^ArAy 



7.13 
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Finally, after setting p *= 1 and Vqo ** 1, for an element 

i s . 



= — = — 7.14 

P AxAy Ax 

This provides the ACp for each element, which can be plotted 
versus chord, for each span-wise section. 

b. Finding the Lift, Drag and Moment 

Lift, drag and moment are also found from the 
vector force on each element. Equation 7.10 showed that the 
force on an element has components normal and tangent to the 
wing. Once more setting p = 1, and Vco = 1, the vector 

force on an element is. 



^F = 



w ~ sin a lATAyf + cos a AT k 



1 . 15 



These forces are oriented relative to the wing coordinate 
system, where i is parallel to the wing and k is 

perpendicular to the wing. Lift and drag are normal and 
tangent, respectively, to the remote velocity (V^o) . The 
wing is at an angle of attack (oc) to that remote velocity. 
Using this information, the program transforms the forces to 
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a coordinate system oriented to the remote velocity. For an 
individual element the lift and drag forces are; 

LL — ElemenLal LiflForce — (cosa)(t^y^V)(cosa) — {iv — sin a) (^y AD {sin a) 

AD = Elemental Drag Force = {cos a) {AyAD (sin a) + (u; — sm a) (AyAD(cosa) 

using a small angle approximation, 7.16 

AL — Elemental Lift F orce — Ay AT — LuaAyAF = Ay AT 
AD = Elemental Drag Force = AyAFa + — a)AyAF = u;AyAF 

The total lift and total drag acting on the wing is a 
summation of the elemental lifts and drags. Using wing 
symmetry ; 



L = Total Lift — 2 ^ AL 

semi — span 

D — Total Drag — 2 ^ AD 

semi —span 



1 . 17 



.The moment about the y axis for an individual 
element, with the accepted sign convention, is; 



AM = Elemental Moment — —{normal wing force) {x position of center of element) 



7.18 



AM = Elemental Moment — ( — AyAF)(r ) 

cen 
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Total moment about the y axis is a summation of 
the elemental moments. Once again, using wing symmetry. 



M — Total Moment — 2 ^ AM 

semi — span 



7.19 



6 . Finding the Coefficients 

Non- dimens ional iz ing the total lift and drag gives, 






L 



-pV^S 
2 “ 






D 



2 ® 



7 . 20 



where S = span x average chord = (b x c), p = 1 and Voo = 1. 

Non - dimens i onal i 2 ing the total moment about the y 
axis gives , 



1 7.21 

-pV cS 

2 ” 

where c is the average chord. 

The -aerodynamic center is important. It is the 

chord-wise position about which the flat wing has zero 
moment. The program finds from. 



QC 



— total moment 
total lift 



1 .22 
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This assumes that total moment is taken about the y axis, or 
x = 0 . The program finds the moment on the elliptically 
loaded wing, and uses the the flat wing as the 
reference axis to find . 

E. TECHNIQUE FOR SPECIFIED LOADING 

1 - Finding Camber and Twist from Loading 

Elliptic loading in a span-wise direction results in 
minimum induced drag, and loading is proportional to 
circulation. At sufficiently high aspect ratio, a flat wing 
with elliptic area distribution generates elliptic loading 
but weighs more than a straight wing. Wings with straight 
leading and trailing edges also cost less to manufacture 
than elliptic wings. So, a span-wise elliptic lift 
distribution is generated by slightly twisting the wing. 
The direct relationship between wing shape (camber and 
twist) and circulation distribution makes this generation 
possible . 
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camber 




Camber and Twist on the Elliptically Loaded Wing 

Figure 7-5 

For a thin wing, elliptic load distribution in the 
chord-wise direction requires an approximately parabolic 
camber shape. The program specifies an elliptic chord-wise 
load distribution and then finds the associated shape. 

2 . Specifying the Lift Coefficient 

Specifying the form and amplitude of the ideal lift 
distribution over the wing determines the form and amplitude 
of the wing shape. A load distribution of elliptic form is 
imposed in this case and the corresponding amplitude is 
fixed by the desired ideal wing lift coefficient, . 

3 . Maximum Desired Lift Coefficient 

The small angle approximation breaks down when large 
lift coefficients are specified. To prevent this 

occurrence, the maximum desired lift coefficient is 
res t r icted to 0.5. 
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4 . Achieving the Desired Lift Coefficient 

The program scales a reference distribution of ACp , 
which is elliptic in form, to achieve the desired total lift 
coefficient for the wing. 

5 . Forces. Moments and Coefficients 

Once the ACp distribution is scaled for the 
elliptically loaded wing, the program finds forces, moments 
and coefficients in the same manner as for the flat wing. 
The only difference is in the coordinate transformation used 
to convert forces on the wing to lift and drag. The 
individual elements on the elliptically loaded wing are no 
longer at a uniform angle to the remote velocity. This 
change in angle must be taken into account when performing 
the coordinate transformation which gives lift and drag 
forces on individual elements. 

The program finds moment coefficient about the 
aerodynamic center, , from. 



Mac Mo Li 

where ^ac is found from the flat wing, c is the average 
chord, C^o pitching moment of the elliptically 

loaded wing about the y axis, and is the lift 

coefficient of the elliptically loaded wing. 



QC 

c 



1 ,23 
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F. 



FLOW CHART FOR MODULES 
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INPUT VARIABLES 



G . 

There are eight input variables. They are: 

1. Aspect Ratio. AR , defined as b^/S, where h is the 

total span and S is the total area, can have any 

positive value less than 20.0. Higher aspect ratios 
should use lifting line theory. 

2. Taper Ratio. LAMDA , defined as where is the 

wing tip chord and C is the wing root chord, can have 
any positive value between 0 and 1. 

3. Leading Edge Sweep. DELTA, defined as the angle in 
degrees that the leading edge of the wing makes with a 
line perpendicular to the remote velocity, can have 
any positive value between 0 and 60. 

4. Angle of Attack. ALPHA, defined as the angle in 

degrees made between the flat wing and the remote 
velocity, can have any positive value between 0 and 
10 . 

5. Number of elements in a section. NX, defined as the 
integer number of elements in a chord of the wing, can 
have any value between 1 and 10. 

6. Number of elements in a semi-span. NY, defined as the 
integer number of elements in a semi-span of the wing, 
can have any value between 1 and 10. 

7. Desired lift coefficient. CLDSRD , the desired lift 
coefficient for the elliptically loaded wing, can have 
any value between 0 and 0.5. 
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H. OUTPUT VARIABLES 



1 . Tabular Data 

When complete, the program writes four tabular data 
files. One file, the SECTION. MAT file, contains data for 
each section in the wing semi-span. Another, the WING. MAT 
file, contains data for the individual grid elements on the 
wing semi-span. The third, the FLAT.DAT file, shows the 
coefficients for the flat wing. The fourth, the CAMBER.DAT 
file, shows the coefficients for the cambered or 
elliptically loaded wing. 

a. SECTION. MAT 

Each row, or line, in the file represents a 
section of the wing semi-span. Each column represents a 
different variable within the section. 



Renofe Velocity 

> 




delta 



y 



A Single section in the Wing Semi-Span 
Figure 7-6 
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The columns are: 



1. Section chord length 

2. Section lift coefficient per unit span for the flat 

wing. The integral of these values from tip to tip is 
the total lift coefficient. 

3. Section drag coefficient per unit span for the flat 

wing. The integral of these values from tip to tip is 
the total drag coefficient. 

4. Section moment coefficient per unit span about the 

aerodynamic center for the flat wing. The integral 

of these values from tip to tip is the total moment 
coefficient about the aerodynamic center. 

5. Section » aerodynamic center relative to x==0 , for 

the flat wing. 

6. Section delta y 

7. Section lift coefficient per unit span for the 

cambered wing. The integral of these values from tip 
to tip is the total lift coefficient, 

8. Section drag coefficient per unit span for the 

cambered wing. The integral of these values from tip 
to tip is the total drag coefficient. 

9. Section moment coefficient per unit span about the 
aerodynamic center for the cambered, or elliptically 
loaded wing. The integral of these values from tip 
to tip is the total moment coefficient about the 
aerodynamic center. 
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10. Section twist in degrees of the cambered wing. 



Positive angle of twist is nose up. 



The 


program 


arranges 


the 


rows 


from 


root to 


tip , with 


the 


root 


section on 


the 


first 


line 


and 


the tip 


section on 


the 


last 


line . 




















b. 


WING 


.MAT 


















Each 


row , 


or 


line 


in the file, 


represents 


an 



individual element in the wing grid. Each column represents 
a different variable for that element. Figure 7-7 shows the 
location of some variables on a representative grid element. 




(x,yb) 

(xp,yp) (force) 
Cxp^yp) (shape) 

Cx,ya) 



A Single Element in the Wing Semi-Span 
Figure 7-7 
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The columns are: 



1. Sequential integer identifier of the element, also 

called IJ in the program. 

2. Integer y position of the element. 1 is at the wing 
root, NY is at the tip. 

3. Integer x position of the element. 1 is at the 
leading edge, NX is at the trailing edge. 

4. X position of center of element. 

5. y position of center of element. 

6. X position of horse shoe vortex, a field point. 

7 . Blank , or zero . 

8. ya position of one corner of horse shoe vortex. 

9. yb position of one corner of horse shoe vortex. 

10. xp trailing edge of the element, a control point used 
in computing wing shape or AT. 

11. xp center of the element, a control point used in 

computing the forces acting on the wing. 

12. yp center of the element, a control point. 

13. Delta x, chord-wise dimension of individual element. 

14. Delta y, span-wise dimension of individual element. 

15. Flat wing incremental circulation (AF). 

16. Flat wing ACp . 

17. Cambered wing incremental circulation (AF), adjusted 

for C^ref. 

18 . Cambered wing ACp 

19. Cambered wing slope (3z/3x). 
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20. Cambered wing height (z). 

21. Cambered wing incremental circulation (AT), if Cj^ref 
were = 1.0. 

22. Down-wash (w) at center of flat wing element due to 
trailing filaments. 

23. Down-wash (w) at center of cambered wing element due 
to trailing filaments. 

Some columns are identical. Separate columns were used for 
each variable during the development phase to promote ease 
in modification. The rows, or lines, of the WING. MAT file 
are arranged as follows. 



Remote Velocity 

> 



N-2 




Grid Element Numbering on the Semi-Span 
Figure 7-8 

The first row is the leading edge element in the root 
section. The next element is immediately behind the leading 
edge and so on. The last row is the trailing edge element in 
the wing-tip section. 
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FLAT.DAT 



c . 



This file contains coefficients and results for 



the flat wing. An example is shown. 



FLAT WING 
CL 




.177721 


CD 


= 




.006516 


CD/CL2 


= 




.206312 


CMAC 


= 




.000000 


XAC 


= 




.279935 


AR 


= 




1.330000 


LAMDA 






.500000 


DELTA 


= 


25.000000 


NX 


= 


8 




NY 


= 


8 




ALPHA 


= 




5.000000 


CLDSRD 


= 




.200000 


ELLIP 


- YES 




CLa, Lift 


Curve Slope 


per/Deg= 




.035544 



d. CAMBER.DAT 



This file contains 
the cambered, or elliptically 
shown . 



coe f f ic ients 
loaded , wing . 



and results for 
An example is 



CAMBERED 


WING 


CL 




.199311 


CD 




.009040 


CD/CL2 




.227557 


CMAC 


= 


- .062732 


AR 


= 


1.330000 


LAMDA 


= 


.500000 


DELTA 


= 


25.000000 


NX 




8 


NY 


= 


. 8 


ALPHA 


= 


5.000000 


CLDSRD 


= 


.200000 


ELLIP 


= 


YES 
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SAMPLE PROBLEM 



I . 



A sample problem will illustrate use of the VORTEX 
program. A wing planform with the following char ac t er is t ic s 
will be analyzed. 



Aspect ratio 
Taper ratio 
Leading Edge Sweep 
Angle of Attack 
Number of span elements 



1.333 
0 . 5 

25 degrees 
5 degrees 
8 



Number of chord elements 8 
Desired Lift Coefficient 0.2 
The planform looks like Figure 7-9. 



Renofe Velocity 

> 



0.5000 -H 




Planform used in the Sample Problem 
Figure 7-9 
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1 . 



Startinp^ the Prograjm 



After turning 
which generally looks 
C :> 

Change the program to 
CD\VORTEX [return] 
the screen should now 



the computer on go to the DOS 
1 ike this . 

the VORTEX directory by typing 
look like. 



p r omp t , 



C\VORTEX:> 

To start the program, type 
VORTEX [return] 

The program will start and the following menu will be 
displayed . 



MAIN PROGRAM MENU 

1. Set Wing Planfo3:m, compute Coefficients and Shape 

2. Set Graphic Display 

3. View Graphic Display 
k . End the Program 



Use this option to set up 
your wing* planform. It 
will also compute all the 
coefficients you need. 
Lift, Drag, Moment, as 
well as the shape, if you 
want an elliptic load 
distribution - 



Main Program Menu 
Figure 7-10 

2 . Main Program Menu 

There are four options in the MAIN PROGRAM MENU. 
All coefficients and results are computed in option 1. 
Whenever the planform is changed, option 1 must be selected 



85 



to recompute the coefficients. The coefficients and results 
are saved to data files after being computed in option 1. 
The data files allow you to use options 2 and 3 as often as 
desired without re-computing the results for the planform. 

Help menus are available. Press FI to toggle the 
help menus on or off. 

3 . Planform and Coefficient Menu 

The program saves planform variables and results 
from the most recent solution and shows them in the PLANFORM 
and COEFFICIENT MENU. The menu looks like Figure 7-11. 



PLANFORM AND COEFFICIENT MENU 

1. Aspect Ratio 

2. Taper Ratio 

3. Leading Edge Sweep Angle (degrees) 

4. Ancle of Attack (degrees) 

5. Number of Elements in a Section 

6. Number of Element.s in a Semi -Span 

7. Compute Wing Shape with Load 

8. E'esired Lift Coefficient 

9. Compute Cofficients for this Planform 

10. Return to Main Menu 



Aspect ratio has values 
between 0 and 10. 



Planform and Coefficient Menu 
Figure 7-11 

Options 1, 2, and 3 change Aspect Ratio, Taper 

Ratio, and Leading edge sweep angle, respectively. 

Option 4 changes the angle of attack. 

Option 5 changes the number of chord-wise grid 
e 1 ements , and 



1 .330000 
. 500000 
25 . 000000 
5.000000 
8 
8 

YES 

. 200000 
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Option 6 changes the number of span-wise grid 
elements in a semi-span. 

Option 7, Compute Wing Shape with Load, is YES or 
NO. Because the time required to compute wing shape is 
extensive, that option may not be desirable for all 
problems . 

Option 8 changes the desired lift coefficient, CLj[ 
of the elliptically loaded wing. 

Option 9 computes the coefficients. The 

coefficients are written to four files, SECTION. MAT, 
WING. MAT, FLAT.DAT, and CAMBER.DAT. Section H of this 

thesis contains a description of the elements in the files. 

Option 10 returns to the Main Program Menu to view 
the results . 
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4 . Graphics Proeram Menu 



Choosing selection 2 of the Main Program Menu 
displays the Graphics Program Menu, An example is shown in 
F igure 7 - 12 . 



GRAPHICS PROGRAM MENU 



1. Type of Wing 

2 . Section Number 

3 . Return to Main Menu 

Spanwise Plots 

4 . d ( CI.i ) /dy 

5. d(CD)/dy 

6. d(CMAC)/dy 

7 . Camhei'ed Wing Twis 

Chordwise Plot 

8. Delta Cp 

9. Wing Shape 



CAMBERED 

5 



Use this option to return 
to the Main Program Menu, 
where you can view the 
graph you have set up 
here . 



Graphics Program Menu 
F igure 7-12 

Option 1 changes the type of wing that will be 
displayed, either flat or cambered. 

Option 2 changes the section of the wing which is 
plotted with options 8 and 9. . 

After setting options 1 and 2, select the desired 
plot from options 4 through 9. The necessary data files 
will be prepared and the program will return to the Main 
Program Menu, where the selected plot can be seen using the 
View Graphic Display option. Examples follow. 
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ft 



FbM WING 

CL = .1693, AOA = 5.00 deg 




FLAT WING d(CL)/dy vs SPAN 
Figure 7-13 





CAMBERED WING 




CLi = .1993, CLi = .20 
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SPAN-WISE COORDINATE Y 



CAMBERED WING d(CL)/dy vs SPAN 
F i gur e 7-14 
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FLAT WING 

CD = .0062. AOA = 5.00 deg 




FLAT WING d(CD)/dy vs SPAN 
Figure 7-15 



CAMBERED WING 

CD = .0090, CLi = .20 




CAMBERED WING d(CD)/dy vs SPAN 
Figure 7-16 
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FLAT WING 

CMAC = .0000. AOA = 5.00 deg 




FLAT WING d(CMAC)/dy vs SPAN 
Figure 7-17 





CAMBERED WING 
CMAC = -.0625. 


CLi = .20 




— U-UUU 

0. 


D O.'l 0.2 


d 
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d 


i 

0.5 
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\-0.036 - 
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< 
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u 








-0.072 - 








-0.090 


SPAN-WISE 


COORDINATE Y 





CAMBERED WING d(CMAC)/dy vs SPAN 
F i gur e 7-18 
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Fb\T WING, SECTION 1 




FLAT WING DELTA CP vs CHORD, Section 1 
Figure 7-19 



CAMBERED WING. SECTION 1 




CAMBERED WING DELTA CP vs CHORD, Section 1 

F igure 7-20 
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PUT WING, SECTION 4 




FLAT WING DELTA CP vs CHORD, Section 4 
F i gure 7-21 



CAMBERED WING, SECTION 4 




CAMBERED WING DELTA CP vs CHORD, Section 4 

F igure 7-22 
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FLAT WING. SECTION 7 




FLAT WING DELTA CP vs CHORD, Section 7 
F i gur e 1-21 



CAMBERED WING, SECTION 7 




CAMBERED WING DELTA CP vs CHORD, Section 7 

Figure 7-24 





FLAT WING, SECTION 8 
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Figure 7-25 




CAMBERED WING, SECTION 8 
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CAMBERED WING DELTA CP vs CHORD, Section 8 

Figure 7-26 
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CAMBERED WING. SECTION 1 
vertical scale exaggerated 




CAMBERED WING SHAPE vs CHORD, Section 1 
F i gur e 7-27 



CAMBERED WING. SECTION 4 
vertical scale exaggerated 




CAMBERED WING SHAPE vs CHORD, Section 4 
Figure 7-28 
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CAMBERED WING. SECTION 7 
vertical scole exoggeroted 




CAMBERED WING SHAPE vs CHORD, Section 7 
Figure 7-29 



CAMBERED WING. SECTION 8 
vertical scale exaggerated 




CAMBERED WING SHAPE vs CHORD, Section 8 
Figure 7-30 
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CAMBERED WING 
CU = .20 




CAMBERED WING TWIST vs SPAN 
Figure 7-31 
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J. METHOD FOR CHANGING LOADING DISTRIBUTION 

This is an optional section and may be skipped without 
loss. Elliptic loading is generated in the SET85 

subroutine. The elements which contain the specified 

loading are WING(IJ,21). If you wish to change the loading 
from elliptic, realize that the grid elements are not 
uniform in size. 

When generating a ACp function, consider the shape. The 
load at the wing tips and trailing edge must go to zero. 
The load at the leading edge should preferably be zero to 
avoid a singularity at that point. 

If you change any of the subroutines, the program must 
be re-compiled with a FORTRAN compiler to include the 
changes into the VLAT85.EXE file. The original version of 
the program was compiled by a MICRO-SOFT 4.2 compiler using 
the following commands. 

FL/FPc VLAT85.FOR 
and 

FL/FPc GRAP85.FOR 

If you have a math co-processor, the FPc option allows the 
program to use the co-processor but does not require one. 

K. METHOD FOR CHANGING THE NUMBER OF GRID ELEMENTS 

This is also an optional section and can be skipped. In 
its original form, the program works with up to 100 grid 
elements in the wing semi-span. The program also has a 
maximum of 10 chord sections. Fewer grid elements or chord 
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sections will run without modification. If you desire more 
than 100 grid elements or 10 chord sections, you must modify 
the program. The changes are not extensive, but require 
that the programs be recompiled, as described in section J 
above . 



Increase grid size in square increments. Modify the 
matrix elements as follows. N is the square dimension and 
corresponds to NX or NY. 



For 


example , if 


you 


want 12 


wing 


sections 


or 12 


grid 


elements 


in each chord 


sec t ion , 


make 


the grid 


12x12 


(144 


elements). In this 


example , N = 


12, and (N*N)+1 


= 145. 


The 


matrices 


that must 


be 


changed 


are 


listed below . 


The 


required 


dimens ions 


are 


shown . 











A( (N*N)+1 , (N*N)+1 ) 
AA( (N*N)+1 , (N*N)+1) 
AS ( (N*N) + 1 , (N^N) + 1) 
WING( (N*N)+1 , 23) 
SECTN(N+1 , 10) 



IOC 






These are the line numbers in the programs which 



contain 


matrix 


variables 


that must 


be changed 










Variable 






Program 


A() 


AA() 


AS() 


WINGO 


SECTN( ) 


VLAT85 


76 


77 


77 


74 


74 


LOAD85 


56 


57 


57 


54 


54 


SHAP85 


54 


55 


55 


52 


52 


GAUSS 


19 


- 


- 


- 


- 


MULTI 


19 


20 


20 


17 


17 


MULT2 


19 


20 


20 


17 


17 


MULT3 


19 


20 


20 


17 


17 


GRAP85 


52 


53 


53 


50 


50 


DNWS85 


- 


- 


- 


23 


23 


SET85 


_ 


_ 


_ 


69 
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Before increasing the grid density, consider the 
tradeoff of memory requirement and solution speed. The 10 X 
10 grid contains 100 elements and 100 unknowns. A 12 x 12 
grid contains 144 elements and 144 unknowns. A 20 x 20 grid 
contains 400 elements and 400 unknowns. More than 4 times 
the number of computations are needed to solve a 400 x 400 
versus a 100 x 100 matrix. You should also increase 
precision if you increase the number of grid elements. All 
elements will be smaller and computer accuracy may be less 
than the distances between points in the grid. 
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APPENDIX A. SOURCE CODE LISTINGS FOR VORTEX PROGRAM 



VLAT85 

^ * 

^ PROGRAM VLAT85 * 

•k k 

* THIS IS THE MAIN DRIVER FOR A NUMBER OF SUBROUTINES THAT ARE 

* LINKED TOGETHER TO FORM A PROGRAM THAT COMPUTES THE FORCES * 

* ON A PLANFORM IN INVICID, INCOMPRESSIBLE FLOW. * 

* THE SUBROUTINES CALLED ARE LISTED HERE AND AT THE END OF THE * 

* PROGRAM AS INCLUDE FILES. * 

k k 

* THIS PROGRAM IS A HIGHLY MODIFIED VERSION OF A PROGRAM IN * 

* 'AN INTRODUCTION TO THEORETICAL AND COMPUTATIONAL AERODYNAMICS' * 

* BY JACK MORAN, WILEY AND SONS, NEW YORK, 1984, FOUND ON PAGE * 

*151. * 

•k -k 

* PROGRAM VARIABLES ARE: * 

* AR : ASPECT RATIO * 

* LAMDA : TAPER RATIO * 

* ADELTA : LEADING EDGE SWEEP ANGLE, IN DEGREES * 

* DELTA : LEADING EDGE SWEEP ANGLE, IN RADIANS * 

* NX : NUMBER OF GRID ELEMENTS IN A SECTION * 

* NY : NUMBER OF GRID ELEMENTS ALONG A SEMI -SPAN * 
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* 


ALPHA 




CLDSRD 


'k 


ELLIP 


k 




k 


PI 


k 


NEQNS 


k 


A ( ) 


k 




k 


AA ( ) 


k 


AS ( ) 


k 




k 




k 


WING( ) 


k 




k 




k 




k 




k 


SECTNO 


k 




k 




k 




k 


IIN 


k 


I OUT 


k 


JOUT 


k 


KOUT 


k 


LOUT 



ANGLE OF ATTACK, IN DEGREES * 
DESIRED LIFT COEFFICIENT WITH ELLIPTIC LOADING * 
CHARACTER FLAG TO COMPUTE THE SHAPE FROM AN * 
ELLIPTIC LOADING * 
IS PI * 
NX*NY, THE TOTAL NUMBER OF EQUATIONS * 
A COEFFICIENT MATRIX, USED FOR THE FLOW TANGENCY * 
GETS TRANSFORMED IN GAUSS * 
THE SAME AS A BUT NOT TRANSFORMED IN THE PROGRAM * 
ANOTHER COEFFICIENT MATRIX, USED FOR THE DOWN-WASH * 
AT THE CENTER OF THE GRID ELEMENT. REQUIRED TO * 
GET THE FORCES ACTING ON EACH GRID ELEMENT. * 



A LARGE ARRAY HOLDING VARIOUS IMPORTANT COORDINATES * 
FOR EACH GRID ELEMENT, AS WELL AS THE RESULTS FOR * 
EACH GRID ELEMENT. A COMPLETE DESCRIPTION OF EACH * 
ELEMENT IN THE ARRAY IS CONTAINED IN THE SUBROUTINE * 



SET85 * 
A LARGE ARRAY HOLDING SOME DIMENSIONS FOR THE * 
SECTIONS, AS WELL AS THE COEFFICIENTS FOR EACH * 
SECTION. A COMPLETE DESCRIPTION OF EACH ELEMENT * 



IN THE ARRAY IS CONTAINED IN THE SUBROUTINE SET85. * 



SET TO 5, FOR KEYBOARD INPUT * 
SET TO 6 FOR SCREEN OUTPUT * 
SET TO 10 FOR WING ARRAY INPUT AND OUTPUT * 
SET TO 11 FOR SECTN ARRAY INPUT AND OUTPUT * 



SET TO 12 FOR LAST PARAMETER INPUT AND OUTPUT 
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* 


MOUT : 


SET TO 13 FOR FLAT PLATE WING COEFFICIENTS 


* 


* 


NOUT : 


SET TO 14 FOR CAMBERED WING COEFFICIENTS 


k 


* 


IMOD : 


A FLAG. IF THE PARAMETERS ARE UNCHANGED FROM THE 


k 


* 




MOST RECENT CALCULATION OF COEFFICIENTS, SET = 0 


k 


* 




IF ANY PARAMETERS ARE CHANGED, SET ' 1 


k 


'k 


ICHOIC : 


MENU INPUT VARIABLE 


k 


* 


ESC : 


THE ESCAPE CHARACTER (27) 


k 




ITYPE : 


INTEGER, 1 IF FLAT, 2 IF ELLIPTIC 


k 


'k 


ISECT : 


INTEGER, SECTION NUMBER TO BE PLOTTED. 


k 


* 






k 


* 


SUBROUTINES 


CALLED ARE: 


k 


•k 


CDAT85 , 


FORTRAN, GENERAL DATA ENTRY, REAL 


k 


k 


CDIT85 , 


FORTRAN, GENERAL DATA ENTRY, INTEGER 


k 


k 


LOAD85 , 


FORTRAN, COMPUTE LOAD DISTRIBUTION FROM A GIVEN 


k 


k 




WING SHAPE. 


k 


k 


SHAP85 , 


FORTRAN, COMPUTE SHAPE FROM A GIVEN LOAD 


k 


k 




DISTRIBUTION. 


k 


k 


DELAY , 


FORTRAN, SLOW DOWN THE PROGRAM SO ERROR MESSAGES 


k 


k 




CAN BE SEEN BEFORE BEING ERASED. 


k 


k 






k 



*********************************************************************** 
PROGRAM VLAT85 
REAL LAMDA 
CHARACTER* 3 ELLIP 
CHARACTER*! ESC 

COMMON PI,B,WING(101,23) , SECTN(11 , 10) , AR, LAMDA, DELTA, NX, NY, 
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* ALPHA, CLDSRD,ELLIP,XAC,CLA 



COMMON/COF/A(101,101) , NEQNS 
COMMON/ CAM/AA( 101, 101 ) , AS(101,101) 

COMMON/FI LS/I IN , lOUT , JOUT , KOUT , LOUT , MOUT , NOUT 

IMOD = 1 

AR = 2.0 

LAMDA - 1.0 

ADELTA= 10.0 

NX = 4 

NY = 4 

ALPHA =5.0 

ELLIP = 'YES' 

CLDSRD= 0.5 

ITYPE = 1 

ISECT = 1 

XAC = 0.0 

PI = 2.0*ACOS(0.0) 

ESC = CHAR(27) 



* 

* SET INPUT/OUTPUT VARIABLES AND OPEN THE NECESSARY FILES * 

* MOUT AND NOUT ARE NOT OPENED UNTIL THE COEFFICIENTS HAVE BEEN * 

* COMPUTED . ^ 

"k k 



kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 






f 



IIN = 5 



lOUT = 6 



JOUT = 10 
KOUT =11 
LOUT =12 
MOUT =13 
NOUT = 14 
OPEN(UNIT = IIN) 

OPEN (UNIT = lOUT) 

OPEN (UNIT = JOUT, FILE= ' WING . MAT ' ) 

OPEN (UNIT = KOUT, FILE=' SECTION. MAT' ) 

OPEN(UNIT = LOUT, FILE=' PARAMS . DAT' ) 

10 FORMAT (I6/3(F6.2,/),2(I3/),F6.2/A/F6.2/I2/I2/F6.2) 

30 FORMAT (23F10.6) 

40 FORMAT (10F10.6) 

50 FORMAT (' THAT IS NOT A VALID OPTION. ENTER AN OPTION BETWEEN 1' 

* ' AND 11') 



60 


FORMAT 


(' 


' ,A,A,F10 


62 


FORMAT 


(' 


' ,A,A,I3) 


64 


FORMAT 


(' 


' ,A,A,A) 



'k'k'k'k'k'kic'k'k'k'k'k'k'k'k’k':kic'k'k'k'k'k'k'k'J('k'k':k'Jc'k'Ji:':k'k'k'k'k'k'k'k'kic'k'k'k'k‘k'k'k'k'k'k'k'^ 



'k 

* READ THE LAST SET OF PARAMETERS, IF UNCHANGED, READ WING AND * 

* SECTN DATA * 

* * 



************************************************************************* 
READ (LOUT, 10, ERR=70, END-70) IMOD,AR,LAMDA,ADELTA,NX,NY,ALPHA, 
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* ELLIP,CLDSRD,ITYPE,ISECT,XAC 



70 DELTA = ADELTA*PI/180 . 

B = AR*(1.0+LAMDA)/2.0 
NEQNS = NX*NY 
IF (IMOD.EQ.O) THEN 

READ (JOUT,30) ( (WING(I ,J) ,J=1 , 23) , 1=1 , NEQNS) 
READ (KOUT.40) ( (SECTN(I , J) , J=1 , 10) , 1=1 ,NY) 
ELSE 
ENDIF 



ic 

* CLEAR THE SCREEN AND WRITE THE MAIN MENU. * 

* THIS USES THE ANSI. SYS CONTROL CODES TO CLEAR THE SCREEN AND * 

* POSITION THE CURSOR. IT ALSO USES FLASHUP WINDOWS FOR THE * 

* MENU PROMPTS. * 



iz'ki^ic'ki^'kic'k'k'k'k'k'kic-kic'kic-k-k'Jc'k'k'k'k'kic’k’k-Jc-kicikic'kic'kic'k'Jcic'k'kic’kic'kift-k'k'k'kic'k'ki^'k^ 



100 



WRITE (lOUT,*) ESC,'[2J' 

WRITE (lOUT,*) ''C=ALL/' 

WRITE (lOUT,*) '*W=LOAD/' 

WRITE (IOUT,60) ESC , ' [ 09 ; 69H' , AR 
WRITE (IOUT.60) ESC , ' [ 10 ; 69H' , LAMDA 
WRITE (I0UT.60) ESC , ' [ 11 ; 69H' , ADELTA 
WRITE (IOUT.60) ESC ,'[ 12 ; 69H' .ALPHA 
WRITE (IOUT.62) ESC , ' [ 13 ; 69H' ,NX 
WRITE (IOUT.62) ESC , ' [ 14 ; 69H' ,NY 
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WRITE (IOUT,64) ESC / [ 15 ; 69H ' , ELLIP 
WRITE (IOUT,60) ESC / [ 16 ; 69H ' , CLDSRD 
NEQNS = NX*NY 

■jV * -jV "Ar ■A: * ■Ar •st ■35: * •Ar * •sV * •A’ -sV •sV ^ •A: •sV ■sV’ * - jV ■A’ * * -Ar - jV -A r * -jV -5^ 

•k k 

* READ THE OPTION. IF OUT OF RANGE, PRINT ERROR MSG, AND START * 

* OVER . * 

* k 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

110 READ(IIN,*,END=110,ERR=110) ICHOIC 

WRITE (lOUT,*) ESC,'[2J' 

WRITE (lOUT,*) ESC, ' [00;00H' 

IF ((ICHOIC. LT.l) .OR. (ICHOIC.GT.il)) THEN 
WRITE (IOUT.50) 

CALL DELAY(6) 

GO TO 100 
ELSE 
END IF 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

k k 

* MAKE ALL CHANGES TO PARAMETERS AS DESIRED. ^ 

* VALUES PASSED TO CDAT AND CDIT ARE LOW, HIGH, VARIABLE, AND * 

* NAME OF THE VARIABLE. IMOD IS RETURNED WITH VALUE 1. * 



*********************************************************************** 



IF (ICHOIC. EQ.l) CALL CDAT85(0 . 0 , 20 . 0 , AR, 
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* 'ASPECT RATIO 



' ,IMOD) 



IF (ICHOIC.EQ.2) CALL CDAT85(0 . 0 , 1 . 0 , LAMDA, 

* 'TAPER RATIO '.IMOD) 

IF (ICHOIC.EQ.3) THEN 

CALL CDAT85(0.0,60.0,ADELTA, 

* 'LEADING EDGE SWEEP ANGLE IN DEGREES '.IMOD) 

DELTA = ADELTA*PI/180. 

ELSE 

ENDIF 

IF (ICH0IC.EQ.4) CALL CDAT85(0 . 0 , 10 . 0 .ALPHA, 

* 'ANGLE OF ATTACK IN DEGREES '.IMOD) 

IF (ICHOIC.EQ.5) CALL CDIT85(0 , 10 .NX, 

* 'NUMBER OF ELEMENTS IN A SECTION '.IMOD) 

IF (ICHOIC.EQ.6) CALL CDIT85(0 , 10 ,NY, 

* 'NUMBER OF SECTIONS IN A SEMI -SPAN '.IMOD) 

IF (ICHOIC.EQ.8) CALL CDAT85 (0 . 0 , 0 . 5 , CLDSRD , 

* 'DESIRED LIFT COEFFICIENT '.IMOD) 

IF (ICHOIC.EQ.7) THEN 

IMOD = 1 

IF(ELLIP.EQ. 'YES' ) THEN 
ELLIP = 'NO ' 

ELSE 

ELLIP = 'YES' 

ENDIF 

ELSE 

ENDIF 
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* RUN THE MAIN PART OF THE PROGRAM. COMPUTE THE LOAD, FROM THE * 

* FLAT PLATE SHAPE, AND THEN THE SHAPE, FROM THE ELLIPTIC LOAD * 

* IF THE ELLIP FLAG IS SET. * 

* * 
*****************************************************■*■*■*■**•*■****'*■******* 



115 IF (ICHOIC.EQ.9) THEN 
IMOD = 0 

WRITE (lOUT,*) ESC,'[2J' 

WRITE (lOUT,*) ESC, ' [00;00H' 
WRITE (lOUT,*) ''C=ALL/' 

WRITE (lOUT,*) '‘W=WAIT/' 

CALL LOADS 5 

IF (ELLIP. EQ. 'YES' ) CALL SHAP85 
ELSE 
ENDIF 



************************************************ ■A’**************'*’****'*'**** 
* * 



* OPTION 10, RETURN TO THE MAIN MENU. DO A NORMAL END. WRITE THE * 

* PARAMETERS AND DATA TO FILES AND CLOSE THEM. * 

* * 



*********************** ********-!!r*******'jSr***'!Sr***')!r***'*-**'***********'****** 

IF (ICHOIC.EQ.IO) THEN 
REWIND (JOUT) 

REWIND (KOUT) 
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REWIND (LOUT) 



WRITE (JOUT, 30) ((WING(I.J) ,J=1,23) ,I=1,NEQNS) 

WRITE (KOUT, 40) ((SECTN(I.J) ,J=1,10) ,1=1, NY) 

WRITE ( LOUT , 10 ) IMOD , AR , LAMDA , ADELTA , NX , NY , ALPHA , ELLI P , 
* CLDSRD,ITYPE,ISECT,XAC 
CLOSE (JOUT) 

CLOSE (KOUT) 

CLOSE (LOUT) 

GOTO 1000 
ELSE 
ENDIF 
GO TO 100 

1000 WRITE (lOUT,*) ESC,'[2J' 

WRITE ( lOUT , *) ESC , ' [ 00 ; OOH ' 

WRITE (I OUT,*) ''C=ALL/' 

WRITE (lOUT,*) '‘W=MAIN/' 

END 



THESE ARE ALL THE SUBROUTINES THAT ARE LINKED, 



'k 

'k 

'k 



-k'k'k'k'k-k'k'k'k-kif^'kic'k'k'k'kic'kic'kic'kic'ki^'kic-k'kiic'k-k'k'kiic'kiici'ci^ic'k'k'kic'kic'k'kic'kic'k'kic'k'k'k'k'k-k'k-k'k'k'k'k 

$ INCLUDE: 'SETS 5. FOR' 

$ INCLUDE : ' DNWS 8 5 . FOR ' 

$ INCLUDE: 'GAUSS. FOR' 

$ INCLUDE: 'CDAT85.FOR' 
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$INCLUDE: 'CDIT85.FOR' 
$ INCLUDE : ' LOAD 8 5 . FOR' 
$ INCLUDE: 'SHAP85.FOR' 
$INCLUDE: 'DELAY. FOR' 

$ INCLUDE; 'MULTI. FOR' 
$INCLUDE: 'MULT2.FOR' 
$INCLUDE: 'MULT3.FOR' 
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LOAD85 



'k'kic'kic'k'k'k'kic'k'k'k'k'k'k'k'kic'k'k'k'k'k'k'k'k'ki^'k'k'k'k'k'k'kic'k'k'ki^'k'k'k'k'k'k'kic'k'k'kic'k'k:^ 



•k 



k 



* SUBROUTINE LOAD85 

k 



k 

k 



^ SUBROUTINE TO DEVELOP A LOADING, OR VORTEX DISTRIBUTION, THAT * 

^ IS ASSOCIATED WITH A FLAT PLATE WING. THE LOADING GENERATED ^ 

* IS RELATED TO THE PRESSURE DIFFERENTIAL, DELTA CP, AND HAS THE ^ 

^ SAME BASIC SHAPE. THE SOLUTION IS GAINED BY USING A MATRIX ^ 

* OF INFLUENCE COEFFICIENTS AND THE WING SHAPE, AND SOLVING FOR * 

* THE VORTEX STRENGTHS. * 

k k 

^ ALL REQUIRED PARAMETERS AND VARIABLES ARE CONTAINED IN COMMON. ^ 

* A DESCRIPTION OF THE COMMON BLOCK VARIABLES IS CONTAINED IN ^ 

* OTHER ROUTINES. * 

k k 



k 

k 

k 

k 

k 

k 

k 

k 

k 

k 



VARIABLES USED IN THE SUBROUTINE ARE: 

ALF : ANGLE OF ATTACK IN RADIANS 
DELTA : THE LEADING EDGE SWEEP ANGLE IN RADIANS 
ADELTA : THE LEADING EDGE SWEEP ANGLE IN DEGREES 
CBAR : THE AVERAGE CHORD 

CL : INTERIM VALUE, USED TO COMPUTE A FINAL VALUE 

CD : INTERIM VALUE, USED TO COMPUTE A FINAL VALUE 

CM : INTERIM VALUE, USED TO COMPUTE A FINAL VALUE 

CXJ : INTERIM VALUE, USED TO COMPUTE A FINAL VALUE 

CZJ : INTERIM VALUE, USED TO COMPUTE A FINAL VALUE 



k 

k 

k 

k 

k 

k 

k 

k 

k 

k 
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* CMJ : INTERIM VALUE, USED TO COMPUTE A FINAL VALUE * 

* XAC : X POSITION OF THE AERODYNAMIC CENTER * 

* CLT : TOTAL LIFT COEFFICIENT OF THE FLAT WING * 

* CDT : TOTAL DRAG COEFFICIENT OF THE FLAT WING * 

* CMT : TOTAL MOMENT COEFFICIENT OF THE FLAT WING * 

* ABOUT THE AERODYNAMIC CENTER * 

* CDOCL2 : CD/CL''2, ALSO CALLED K IN THE DRAG POLAR * 

* * 



* SUBROUTINES CALLED ARE: * 

* SET85 : FORTRAN, USED TO GENERATE THE WING AND SECTION * 

* ARRAY CONTAINING PLANFORM COORDINATES. * 

* DNWS85 : FORTRAN, USED TO COMPUTE THE DOWNWASH GENERATED * 

* AT A POINT BY A HORSE-SHOE VORTEX AND IT'S MIRROR * 

* IMAGE ON THE OPPOSITE SEMI -SPAN. * 



* GAUSS : FORTRAN, USED TO SOLVE THE INFLUENCE COEFF MATRIX * 



' * AND SLOPE VECTOR, FOR THE VORTEX STRENGTH. USES * 

j * SIMPLE GAUSS ELIMINATION. * 

* MULT2 : FORTRAN, USED TO COMPUTE THE DOWNWASH GENERATED * 

I 

* ON EACH VORTEX ELEMENT, NECESSARY TO FIND THE * 

* FORCES ACTING ON THE WING. * 

-k k 



kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

SUBROUTINE LOAD85 
REAL LAMDA 
CHARACTER* 3, ELLIP 

COMMON PI,B,WING(101,23) , SECTN(11 , 10) ,AR, LAMDA , DELTA, NX, NY, 
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* ALPHA, CLDSRD,ELLIP,XAC,CLA 



COMMON/COF/A(101,101) , NEQNS 
COMMON/CAM/AA(101,101) , AS(lOl.lOl) 

COMMON/FILS /IIN , lOUT , JOUT , KOUT , LOUT , MOUT , NOUT 
60 FORMAT ( ' FLAT WING', 

* /' CL - ',F12.6/' CD = ',F12.6/' CD/CL2 = ' , F12 . 6 

* /' CMAC = ',F12.6/' XAC = '.F12.6/' AR = ' , F12 . 6 



* /' LAMDA = ',F12.6/' DELTA = ',F12.6/' NX = ',13 

* /'NY = ',13 /' ALPHA = ',F12.6/' CLDSRD = '.F12.6 



* /' ELLIP = ' ,A3 /' CLa, Lift Curve Slope' 

* /' per/Deg= ',F12.6) 

•st -sSr -sSr -it -sSr -A- -sV -5^ -sV "sSr ^ ’A’ -jV - jSr ^ -jV -jV -A" - j?: -sSr ^ ‘sV -jit * ^ ^ 

* ^ 

* COMPUTE THE ANGLE OF ATTACK IN RADIANS. 

•j5r ^ 

ALF = ALPHA*PI/180.0 
SINALF = SIN(ALF) 

COSALF = COS (ALF) 

A- A: A- A: -A -A -A -A -A "A A” "A -A "A A” "A A” -A -A -A -A A- -A "A -A A” •A -A ^ A” A” -A -A A- -A A” •A "A A: A” •A A” ^ •A A- 
A A 

* CALL THE SETUP ROUTINE. USED TO BUILD THE COORDINATES OF THE * 

* WING. THEY ARE STORED IN THE ARRAYS WING AND SECTN. * 

* * 
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CALL SETS 5 

i!c'k'k'k'k'k'k’k'kik'k'k’kik'kik'k-k'k'kic'k'k'k'k'k'k'k'k'k'k'k'k'kic'k'k'k'k'k':k'k'k'k':k':k'kik'k'k'k'k'k'k 



* * 

* BUILD THE THREE COEFFICIENT MATRICES. A(), AA(), AND AS() * 

* A() AND AA() ARE THE SAME, AND ARE FOR EVALUATING FLOW * 

* TANGENCY AT THE EDGE OF THE GRID ELEMENTS. * 

* AS() IS FOR EVALUATING DOWN-WASH GENERATED ON THE VORTEX * 

* ELEMENT ITSELF. * 

* * 



************************************************************************************* 
DO 130 1=1, NY 
DO 120 J=1,NX 

IJ = (I-1)*NX+J 

*********************************************************************** 
* * 

* AUGMENT THE A() MATRIX WITH THE RIGHT HAND SIDE, WING SHAPE. * 

•jV * 

'^k'k'k'k'k'k'k'k'k'k'kic'k'k'k'k'k'kic'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k-k'k'kik-k'^ 

A(IJ,NEQNS+1) = SINALF 

* A(IJ,NEQNS+1) = ALF 
DO 110 K=1,NY 

DO 100 L=1,NX 

KL = (K-1)*NX+L 

CALL DNWS85(IJ,KL,A(KL,IJ) , 'CONT') 
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AA(KL,IJ) = A(KL,IJ) 



CALL DNWS85(IJ,KL,AS(KL,IJ) , 'SELF') 

100 CONTINUE 

110 CONTINUE 

120 CONTINUE 

130 CONTINUE 

*********************************************************************** 
* * 

* COMPUTE THE AVERAGE CHORD * 

* * 
*********************************************************************** 

CBAR = (SECTN(l,l)+SECTN(NY,l))/2.0 

* ”5^ :jSr "st ^ "st 'it ”3^ "st 'jSt'sV -st •sSr 'sSr’st -st "jS: "st 'jSr’jSr ■sV’jSr ^ -jS: ^ "it "it ^ ^ ^ 



■ 35 : * 

* CALL THE GAUSS ROUTINE, TO SOLVE FOR THE LOAD VECTOR * 

* THE RESULT IS RETURNED IN PLACE OF THE ORIGINAL LOAD VECTOR * 

* * 



*********************************************************************** 
CALL GAUSS (1) 

*********************************************************************** 



* * 

* STORE THE RESULT IN THE WING() MATRIX * 

* CONVERT THE CIRCULATION, OR VORTEX STRENGTH IN ELEMENT 15 TO * 

* DELTA CP. IN ELEMENT 16 * 

* * 



*********************************************************************** 



1 ^ 



DO 210 I = 1,NY 



DO 200 J = 1,NX 

IJ = (I-1)*NX+J 
WING(IJ,15) = A(IJ,NEQNS+1) 

WING(IJ,16) = 2.0*WING(IJ,15)/WING(IJ,13) 

200 CONTINUE 

210 CONTINUE 

ic 'fc 

* COMPUTE THE DOWNWASH ON EACH VORTEX ELEMENT, NECESSARY TO * 

* FIND THE FORCES ACTING ON THE WING. * 

* * 
•k'k-kif'ic-k-k-k-k-k-k-k-k-k'k-k-k-k'kicie'ic'k-k-k-k-k'ir'k-k-lrirk-k-k'irk-ie^'k'k-k-k'ir'k-k'k-k-k-k-k'k-k-k-k-k-k-k-k-k-k-k-k-k-lck-k-k-frjdc 

CALL MULT 2 

* ^ 

* INITIALIZE VALUES AND COMPUTE THE FORCES ACTING ON THE WING * 

* SECTIONS AND ALSO THE TOTAL WING. * 

* * 
*********************************************************************** 

CD = 0.0 
CL - 0.0 
CM = 0.0 
DO 330 1=1, NY 
CXJ = 0.0 
CZJ = 0.0 
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CMJ =0.0 



DO 320 J=1,NX 

IJ = (I-1)*NX+J 

CZJ = WING(IJ,15)*2.0*(1. -WING(IJ,22)*SINALF)/(B*CBAR) 

CXJ = WING(IJ,15)*WING(IJ,22)*2.0*COSALF/(B*CBAR) 

CMJ = -WING(IJ,15)*WING(IJ,A)*2.0* 

* (1. -WING(IJ,22)*SINALF)/(B*CBAR*CBAR) 

SECTN(I,2) = SECTN(I,2)+CZJ 
SECTN(I,3) = SECTN(I,3)+CXJ 
SECTN(I,4) = SECTN(I,4)+CMJ 
320 CONTINUE 

CL = CL+SECTN(I,2)*SECTN(I,6) 

CD = CD+SECTN(I,3)*SECTN(I,6) 

CM = CM+SECTN(I,4)*SECTN(I,6) 

SECTN(I,5) = -SECTN(I,4)/SECTN(I,2) 

330 CONTINUE 

* * -jt * -jt -st -jV * ■st "5^ * * * * -3^ ^ * -jt "jt -jV ^ -A- ^ 



■A- -A- 

^ NOW COMPUTE THE MOMENTS ABOUT THE AERODYNAMIC CENTER * 

■k k 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

XAC = -CM*CBAR/CL 
DO 335 I = 1,NY 

SECTN(I,4) = SECTN(I,4) + SECTN(I , 2)*XAC/CBAR 
335 CONTINUE 



*********************************************************************** 
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COMPUTE THE FINAL TOTAL WING COEFFICIENTS AND WRITE TO THE 






•)V 



'k 



k 



FILE. * 

k k 



kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

CLT = 2.0*CL 
CDT = 2.0*CD 

CMT = 2.0*CM + (CLT*XAC/CBAR) 

CLA - CLT/ALPHA 

CDOCL2 = CDT/CLT**2 

ADELTA = DELTA * 180. /PI 

OPEN (UNIT = MOUT, FILE ='FLAT.DAT') 

WRITE (MOUT , 60 ) CLT , CDT , CDOCL2 , CMT , XAC , AR , LAMDA , ADELTA , NX , NY , 

* ALPHA, CLDSRD, ELLIP, CLA 
CLOSE (MOUT) 

340 RETURN 
END 
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SHAP85 



'k'k'k'k'k'k'k'kiciiic'k'kicicic^k'k^k'ki^'k'k'kif'k'kicic'k'k'k'ki^'ki^'ki^ici^'k'ki^i^'kicici^'kic'kic'k'k'k'k'k'k'k'k'k'k'k'k'kic'k'k'k'k'k 
k k 

* SUBROUTINE SHAP85 * 

•k k 

* SUBROUTINE TO DEVELOP A WING SHAPE THAT IS ASSOCIATED WITH ^ 

* A SPECIFIED LOADING FUNCTION. THE LOADING FUNCTION IS USUALLY * 

^ ELLIPTIC, IN THE SPAN AND CHORD DIRECTION. THE SHAPE WILL * 

* USUALLY RESEMBLE A PARABOLIC CAMBER, IN CHORD, AND HAVE SOME ^ 

* DEGREE OF TWIST, ALONG THE SPAN. * 

k k 

* ALL REQUIRED PARAMETERS AND VARIABLES ARE CONTAINED IN COMMON. * 

* A DESCRIPTION OF THE COMMON BLOCK VARIABLES IS CONTAINED IN * 

* OTHER ROUTINES. * 

* * 

* VARIABLES USED IN THE SUBROUTINE ARE: * 

* CBAR : THE AVERAGE CHORD * 

* CLILST : CLREF FROM THE PREVIOUS ITERATION * 

* CLREF : A SCALER, USED TO INCREASE OR DECREASE THE AMPL * 

* OF THE SHAPE, IN ORDER TO ACHIEVE THE DESIRED * 

* LIFT COEFFICIENT. * 

* MOD : A COUNTER USED TO LIMIT THE NUMBER OF ITERATIONS * 

* THAT THE PROGRAM CAN USE IN AN ATTEMPT TO GENERATE * 

* A DESIRED LIFT COEFFICIENT. * 

* H : A TEMPORARY VARIABLE USED WHILE INTEGRATING THE * 

* HEIGHT OF THE WING ALONG A CHORD. * 
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COSCAM : 


THE COSINE OF THE ANGLE MADE BY AN INDIVIDUAL GRID 


k 


'k 




ELEMENT AND THE REMOTE VELOCITY. 


k 


'k 


SINCAM : 


THE SIN OF THE ANGLE MADE BY AN INDIVIDUAL GRID 


k 


'k 




ELEMENT AND THE REMOTE VELOCITY. 


k 


'k 


CCL : 


INTERIM VALUE, USED TO COMPUTE A FINAL VALUE 


k 


ic 


CCD : 


INTERIM VALUE, USED TO COMPUTE A FINAL VALUE 


k 


'k 


CCM : 


INTERIM VALUE, USED TO COMPUTE A FINAL VALUE 


k 


'k 


CXCP : 


X POSITION OF THE CENTER OF PRESSURE 


k 


'k 


CCLT : 


TOTAL LIFT COEFFICIENT OF THE CAMBERED WING. 


k 


k 


CCDT : 


TOTAL DRAG COEFFICIENT OF THE CAMBERED WING. 


k 


k 


CCMT : 


TOTAL MOMENT COEFFICIENT OF THE CAMBERED WING. 


k 


k 




ABOUT THE AERODYNAMIC CENTER 


k 


k 


CCDOCL : 


CD/CL^2, ALSO CALLED K IN THE DRAG POLAR 


k 


k 






k 


k 


SUBROUTINES 


CALLED ARE: 


k 


k 


MULTI : 


FORTRAN, USED TO COMPUTE THE SHAPE OF THE WING 


k 


k 




FROM THE LOADING VECTOR AND THE INFLUENCE COEFF 


k 


k 




MATRIX 


k 


k 


MULT 3 : 


FORTRAN, USED TO COMPUTE THE DOWNWASH GENERATED 


k 


k 




ON EACH VORTEX ELEMENT, NECESSARY TO FIND THE 


k 


k 




FORCES ACTING ON THE WING. 


k 


k 






k 



**********-st*****************************************-jt***-jlr •*■******•*■****■*•■* 

SUBROUTINE SHAP85 
REAL LAMDA 
CHARACTER*3, ELLIP 
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COMMON PI,B,WING(101,23) , SECTN(11 , 10) , AR.LAMDA , DELTA ,NX, NY, 
* ALPHA, CLDSRD,ELLIP,XAC,CLA 
COMMON/COF/A(101 , 101) , NEQNS 
C0MM0N/CAM/AA(101,101) , AS(101,101) 

COMMON/FILS/I IN , lOUT , JOUT , KOUT , LOUT ,MOUT , NOUT 
CLILST =0.0 
CLREF =0.05 
MOD = 0 



66 FORMAT (' UNABLE TO REACH THAT CLDSRD') 



70 


FORMAT 


■ ( ' ( 


CAMBERED WING', 












•A- /' 


CL 


= ' ,F12.6/' CD 


= ' ,F12.6/' 


CD/CL2 = 


' ,F12.6 






^ /' 


CMAC 


= ' ,F12.6/ 


f 


AR 


' ,F12.6 






'k /' 


LAMDA 


= ',F12.6/' DELTA 


- ' ,F12.6/' 


NX 


' ,13 






•A- /' 


NY 


= ' ,13 /' ALPHA 


= ' ,F12.6/' 


CLDSRD - 


' ,F12.6 






•A: /' 


ELLIP 


= '.A3) 












'k 














'k 


•A* 


COMPUTE THE 


AVERAGE CHORD 








•k 


•A- 














'k 



'jV -jt -jS: "jt -jt -jt -A- "sV *5^ ^ "it - jV ' 5^ -jt ”5 ^ -jV 'it “it ^ -5V *5^ -jt -Ar -jt -A- -it -jV “it -jV ”5^ -jt ^ 

CBAR = (SECTN(l,l)+SECTN(NY,l))/2.0 

-A: -A: -A: "A: -A: -A: -A: -A: -Ar -A: "A: ^ -A: -A: -Ar -A: -A: -A: -A: -A: -Ar -A: -A: -A: -A: -Ar -A: Ar -A: "A: -jV Ar -A: -A: -A: -A: -A: -A: -A: 

* * 

* USE THE SCALER, CLREF, TO PUT AN ELLIPTIC LOAD DISTRIBUTION 

* IN THE VECTOR (IJ,17). AT THE SAME TIME COMPUTE THE DELTA CP * 

* IN THE VECTOR (IJ,18). * 



1 



* * 
*********************************************************************** 
390 DO 400 IJ = l.NEQNS 

WING(IJ,17) = WING(IJ ,21)*CLREF 
WING(1J,18) = 2.0*WING(IJ,17)/WING(IJ,13) 

400 CONTINUE 

* * 

* CALL THE SUBROUTINE MULTI WHICH WILL * 

* MULTIPLY THE LOAD VECTOR. (IJ.17) BY THE COEFFICIENT MATRIX, * 

* AA( ), TO GET THE SHAPE VECTOR, (IJ,19) * 

* * 
*********************************************************************** 

CALL MULTI 

*********************************************************************** 

* * 

* CHECK THE RESULTING SHAPE VECTOR (IJ,19), TO SEE IF THE ANGLE * 

* THAT WAS REQUIRED HAS A SIN GREATER THAN 1, THAT IS SIMPLY * 

* CHECKING THE SHAPE VECTOR, TO SEE IF IT IS GREATER THAN 1. * 

* IF IT IS, THE SCALER MULTIPLIER, CLREF, IS DECREASED. AT THE * 

* SAME TIME, MOD, THE ITERATION LIMITER, IS INCREMENTED. IF * 

* MOD IS GREATER THAN 5, THAT MEANS WE OVERSHOT 5 TIMES, AND * 

* PROBABLY WILL NEVER GET TO THE DESIRED LIFT COEFFICIENT. * 

* IN THAT CASE, COMPUTE THE LAST VALID SHAPE, AND PROVIDE THOSE * 

* DATA . * 

* * 
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^:Ar’55r^Vt*Vt:)S:'55r^S:^^^'55r^^-}S::A:'55:^t^'55::A:'35r^'3S:^*'55r*^'35:*'55:^VSr'3S:^VIr-A’*^’:S:’35:VS:*^*^^'5t 



DO 410 IJ = l.NEQNS 

IF (ABS(WING(IJ,19)).GT.1.0) THEN 
CLREF = CLREF*.9 
MOD = MOD + 1 
IF (M0D.GT.5) THEN 
WRITE (IOUT,66) 

CALL DELAY(6) 

GO TO 540 
ELSE 

GO TO 390 
ENDIF 

ELSE 
ENDIF 
410 CONTINUE 

**********************************')S:*******5S:***************')!r***')!r ■)*:*■*•')!:**** 



* * 

* BEGIN TO COMPUTE THE FORCE AND MOMENT COEFFICIENTS ON THE * 

* CAMBERED WING. MULT3 WILL MULTIPLY THE SPECIFIED LOADING 

* BY THE COEFFICIENT MATRIX, AS, TO GET THE DOWNWASH AT THE * 

* CENTER OF EACH CAMBERED ELEMENT. INTERMEDIATE VARIABLES ARE * 

* INITIALIZED. * 

* COEFFICIENTS FOR SECTIONS AND THE WHOLE WING ARE COMPUTED. * 

* * 



'k'k'k'k^'k'k':k'k^ik'k'^^ik'k'^'kik'k'k'k'k':k'k'k'k'k'k'k'k'k'k'k'k'k'ki(:ic'k'k'k'k'k'kiz'^ 

CALL MULT3 
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CCD =0.0 



520 



CCL =0.0 
CCM =0.0 
DO 530 1=1, NY 
CCXJ - 0.0 
CCZJ =0.0 
CCMJ =0.0 
SECTN(I,7) = 0.0 
SECTN(I,8) =0.0 
SECTN(I,9) = 0.0 
H = 0.0 

IJ = (I-1)*NX+1 
DO 520 J=1,NX 

IJ - (I-1)*NX+J 
H = H-WING(IJ,19)*WING(IJ,13) 

WING(IJ,20) = H 

COSCAM = COS(ASIN(WING(IJ,19))) 

SINCAM = WING (IJ, 19) 

CCZJ = WING(IJ,17)*2.0*(1. -WING(IJ,23)*SINCAM)/(B*CBAR) 
CCXJ = WING(IJ,17)*WING(IJ,23)*2.0*COSCAM/(B*CBAR) 

CCMJ = -WING(IJ,17)*WING(IJ,4)*2.0* 

* (1. -WING(IJ,23)*SINCAM)/(B*CBAR*CBAR) 

SECTN(I,7) = SECTN(I,7)+CCZJ 
SECTN(I,8) = SECTN(I,8)+CCXJ 
SECTN(I,9) = SECTN(I,9)+CCMJ 
CONTINUE 
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CCL = CCL+SECTN(I,7)*SECTN(I,6) 

CCD = CCD+SECTN(I,8)*SECTN(I,6) 

CCM = CCM+SECTN(I,9)*SECTN(I,6) 

******************************************************************************** 



* * 

^ COMPUTE THE SECTION TWIST. THE CAMBERED WING CL * 

* WILL BE THE SAME AS CLDSRD . WITHIN A VERY SMALL MARGIN. ^ 

'k k 



kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

SECTN(I,10) = -WING(I*NX,20)*180./(SECTN(I,1)*PI) 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 
k k 

* COMPUTE THE MOMENT COEFFICIENT ABOUT THE AERODYNAMIC CENTER. * 

k k 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

SECTN(I,9) = SECTN(I,9) + (SECTN(I , 7)*XAC/CBAR) 

530 CONTINUE 

CXCP = -CCM/CCL 
CCLT = 2.0*CCL 
CCDT = 2.0*CCD 

CCMT = 2.0*CCM + (CCLT*XAC/CBAR) 

CCDOCL = CCDT/CCLT**2 
CLILST = CLREF 

IF (ABS(CLDSRD-CCLT) .LT. (CLDSRD*0.01)) GO TO 540 
CLREF = CLREF*CLDSRD/CCLT 
GO TO 390 
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'k ^ 

* THE DESIRED LIFT COEFFICIENT HAS BEEN ACHIEVED, WITHIN A * 

* A TOLERANCE OF +/- 1.0%. THE CORRECT VORTEX STRENGTH AND * 

* DELTA CP ARE PUT INTO THE WING ARRAY, AND THE FINAL CAMBERED * 

* WING VALUED OF LIFT, DRAG, K, MOMENT, AND CENTER OF PRESSURE * 

* ARE PRINTED. * 

* * 



540 CONTINUE 

DO 550 IJ = 1,NEQNS 

WING(IJ,17) = WING(IJ,21)*CLILST 
WING(IJ,18) = 2.0*WING(IJ,17)/WING(IJ,13) 

550 CONTINUE 

*********************************************************************** 
* * 

* WRITE THE TOTAL WING COEFFICIENTS TO THE SCREEN * 

* * 
*********************************************************************** 

OPEN (UNIT = MOUT, FILE = 'CAMBER.DAT') 

ADELTA = DELTA *180. /PI 

WRITE (MOUT, 70) CCLT,CCDT,CCDOCL,CCMT, AR,LAMDA, ADELTA, NX, NY, 

* ALPHA, CLDSRD, ELLIP 
CLOSE (MOUT) 

RETURN 

END 
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SET85 



'k'k:k'k:k’^:k:k'^^'k:k:k:k:k:k^'^^'k:kic'k^'k:k'k'k^'kic:k'k'k'k:k:k'k^'k:k:k'k'^'k^'^:k:k'k:k:k:k'^:k^'k'k'k'k-k'k'k:k'k:k:k'k'k'k'k'k 



'k 

'k 



'k 

'k 



SUBROUTINE SET85 

This program is written for the setup of a straighthorse shoe 
vortex over tapered and swept wing. It uses cosine spacing 
for both the span and chord elements. This concentrates 
elements near the edges of the wing, where they are most 
important . 



* 



•k 



•k 

k 



k 

k 



Column elements in the WING array are: 

1 sequential position, also called IJ 

2 integer y position 

3 integer x position 

4 X center of element 

5 y center of element 

6 X position of horse shoe vortex 

7 blank 

8 ya position of one corner of horse shoe vortex 

9 yb position of one corner of horse shoe vortex 

10 xp position of flow tangency point, 3/4 chord 

11 xp position of down wash point, 1/4 chord 

12 yp position of flow tangency or down wash pt, center of elem * 

13 del X, chord wise dimension * 

14 del y, span wise dimension 



r 



p 

< 





15 


flat plate vorticity 






•A- 


•A* 


16 


flat plate delta cp 






k 


■k 


17 


cambered voricity, adjusted for CLi 






k 


•A- 


18 


cambered delta cp 






k 


k 


19 


cambered slope 






k 


k 


20 


cambered height 






k 


k 


21 


cambered vorticity, if CLREF were == 1.0 






k 


k 


22 


downwash at center of flat plate element 






k 


k 


23 


downwash at center of cambered element 






k 


k 










k 


k 


Column elements in the SECTN array are: 






k 


k 


1 


chord 






k 


k 


2 


flat wing lift force per unit span 






k 


k 


3 


flat wing drag force per unit span 






k 


k 


4 


flat wing moment about the aerodynamic center 




k 


k 


5 


flat wing XCP, center of pressure relative 


to x= 


=0 


k 


k 


6 


section delta y 






k 


k 


7 


cambered wing lift force per unit span 






k 


k 


8 


cambered wing drag force per unit span 






k 


k 


9 


cambered wing moment about the aerodynamic 


center 


k 


k 


10 


cambered wing twist. 






k 


k 










k 


k 


VARIABLES USED ARE: 






k 


k 




B : WING SPAN 






k 


k 




DTHE : AN ANGULAR ELEMENT WIDTH, IN RADIANS 


SPAN 


k 


k 




DPS I : AN ANGULAR ELEMENT DEPTH, IN RADIANS, 


CHORD 


k 
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* 


THEA 


• 


INTERMEDIATE VARIABLE 


'k 


'k 


THEA 




INTERMEDIATE VARIABLE 


* 


* 


PSIA 




INTERMEDIATE VARIABLE 




* 


PSIB 




INTERMEDIATE VARIABLE 


'k 


* 


ELE 




INTERMEDIATE VARIABLE 


•k 


* 


ELT 




INTERMEDIATE VARIABLE 


'k 


* 


PSI 




CHORDWISE ANGULAR COORD FOR THE CENTER OF AN 


* 


* 






ELEMENT. USED IN GENERATING ELLIPTIC LOAD. 


* 


* 


THETA 




S PANWISE ANGULAR COORD FOR THE CENTER OF AN 


•k 


* 






ELEMENT. USED IN GENERATING ELLIPTIC LOAD. 


'k 


* 








'k 


* 


FUNCTIONS 


USED ARE: 




* 


BETA 


: 


COMPUTES THE LENGTH OF A LOCAL SECTION CHORD. 


'k 


•jV 








* 


* 


ALL VARIABLES AND PARAMETERS ARE PASSED IN COMMON BLOCKS 


★ 


* 


THE ONLY VARIABLES CHANGED IN THIS ROUTINE ARE: 


'k 


* 


B 






* 


* 


WING() 








* 


SECTNO 






'k 


* 








•k 



'A' •A’ •sV -sV -sV 'sV -sV •A’ •5V •A’ ^ •sV "sV "st -sV ■sV "it -sV * -it -sV * -jV •:5r - jV * * -sV 

SUBROUTINE SET85 

COMMON PI,B,WING(101,23) , SECTN(11 , 10) ,AR ,LAMDA,DELTA,NX,NY , 

* ALPHA, CLDSRD,ELLIP,XAC,CLA 
REAL LAMDA 

*************************************************************** ■*■**■* ********* 
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•k ^ 

* DEFINE THE FUNCTION BETA, WHICH COMPUTES THE LENGTH OF A LOCAL ^ 

* CHORD . 

'3^' * 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk'k'kkk'kkkkkkk'k 

BETA(Y,LAMDA,AR) = 1 . 0- (Y*4 . 0*(1 . 0-LAMDA)/(AR*(l . 0+LAMDA) ) ) 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk'k'k'kk'k'k'kk 



k k 

* DEFINE THE SIZE OF SPAN AND CHORD WISE GRID ELEMENTS. * 

* DETERMINE THE TOTAL WING SPAN. * 

* * 



•k-k-k-k-k-k-k-k-k-k-k-ir/rk-irk-k-kirk-kirk-k-kirk-k-k-kirk-k-kirk-k-kirk-kicfrk-k-k-k-k-kirk-k-kirk-k-kirk'kick-k-k-k-k-kick-ic-k 

* DTHE = PI*0. 5/FLOAT (NY) 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 



k k 

* THIS STATEMENT SETS THE GRID UP SO THAT THE OUTER ELEMENT * 

* IS IGNORED. * 

k k 



kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

DTHE = PI*0. 5/FLOAT (NY+1) 

DPS I = PI/FLOAT (NX) 

B = AR*(1.0+LAMDA)/2.0 
DO 110 I = 1,NY 
DO 100 J = 1,NX 

*************************** •A’*********************************************** 
* * 
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COMPUTE ALL REQUIRED COORDINATES FOR THE WING() ARRAY. 



* 



•A- 



•A- -A- 

IJ = (I-1)*NX+J 
WING(IJ,1) = FLOAT(IJ) 

WING(IJ,2) = FLOAT(I) 

WING(IJ,3) = FLOAT(J) 

THEA = (FLOAT(I-1)*DTHE)+(PI/2.0) 

THEB = (FLOAT(I)*DTHE)+(PI/2.0) 

PSIA = (FLOAT(J-l)*DPSI) 

PSIB = (FLOAT(J)*DPSI) 

WING(IJ,8) = -B*GOS(THEA)/2.0 
WING(IJ,9) = -B*COS(THEB)/2.0 
WING(IJ,12) = -B*(COS(THEA)+COS(THEB))/4.0 
ELE = (1.0 - COS(PSIA))* 

* BETA(WING(IJ,12) ,LAMDA,AR)/2.0 
ELT - (1.0 - COS(PSIB))* 

* BETA(WING(IJ ,12) ,LAMDA,AR)/2.0 

* THESE TWO STATEMENTS SET THE POINTS AT 1/4 AND 3/4 CHORD 

WING(IJ.IO) = WING(IJ,12)*TAN(DELTA)+(ELE+3*ELT)/4.0 
WING(IJ,11) = WING(IJ ,12)*TAN(DELTA)+(3*ELE+ELT)/4.0 
WING(IJ,4) = WING(IJ,12)*TAN(DELTA)+(ELE+ELT)/2.0 
WING(IJ,6) = WING(IJ,11) 

WING(IJ,5) = WING(IJ,12) 

WING (IJ, 13) = ELT- ELE 
WING(IJ,14) = WING(IJ,9)-WING(IJ,8) 
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PSI = ACOS(1.0-(2.0*(WING(IJ,11)-(WING(IJ,5)* 

* TAN(DELTA)))/BETA(WING(IJ,5) ,LAMDA,AR))) 

THETA = ACOS(-2.0*WING(IJ,5)/B) 

WING(IJ,21) = 4.0*(1.0+LAMDA)*SIN(THETA)*SIN(PSI)* 

* WING(IJ,13)/(PI**2*BETA(WING(IJ,5) ,LAMDA,AR)) 

WING(IJ,7) = 0.0 

WING(IJ,15) = 0.0 
WING (IJ, 16) =0.0 
WING(IJ,17) = 0.0 
WING(IJ,18) = 0.0 
WING(IJ,19) = 0.0 
WING(IJ,20) = 0.0 
WING(IJ,22) = 0.0 
WING(IJ,23) = 0.0 
100 CONTINUE 

* * 

* COMPUTE ALL REQUIRED COORDINATES FOR THE SECTN() ARRAY, 

•k k 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk'k 

SECTN(I.l) = BETA(WING(IJ,12) ,LAMDA,AR) 

SECTN(I,2) = 0.0 
SECTN(I,3) =0.0 
SECTN(I,4) = 0.0 
SECTN(I,5) = 0.0 
SECTN(I,6) = WING(IJ,14) 
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SECTN(I,7) =0.0 



SECTN(I,8) = 0.0 
SECTN(I,9) = 0.0 
SECTN(I.IO) =0.0 
110 CONTINUE 
RETURN 
END 



r 






DNWS85 



******************************************************************************** 



'k 






k 


k 


SUBROUTINE 


DNWS85(IJ,KL,W,IND) 


k 


k 






k 


k 


COMPUTE DOWNWASH ON GRID ELEMENT NUMBER KL DUE TO 


k 


k 


VORTICES ON GRID ELEMENT IJ, AND IT'S MIRROR IMAGE 


k 


k 






k 


k 


VARIABLES 


USED ARE: 


k 


k 


IJ 


: VORTEX GRID ELEMENT NUMBER, INPUT 


k 


k 


KL 


: DOWNWASH GRID ELEMENT NUMBER, INPUT 


k 


k 


W 


: THE DOWN WASH, RETURNED TO THE CALLING PROGRAM 


k 


k 




OUTPUT 




k 


IND 


: A FLAG TO DETERMINE WHETHER TO USE TRAILING 


k 


k 




EDGE OF THE GRID ELEMENT, OR THE CENTER OF 


k 


k 




THE ELEMENT, INPUT 


k 


k 






k 


k 


FUNCTIONS 


USED ARE: 


k 


k 


WHVl 


: FORTRAN, COMPUTES THE DOWNWASH DUE TO ONE CORNER 


k 


k 




OF THE HORSESHOE VORTEX 


k 


k 


WHV2 


: FORTRAN, COMPUTES THE DOWNWASH DUE TO JUST THE 


k 


k 




TRAILING FILAMENTS OF THE HORSESHOE VORTEX 


k 


k 






k 



****************************************************************************** 
SUBROUTINE DNWS85 (IJ ,KL,W, IND) 
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COMMON PI,B,WING(101,23) , SECTN(11 , 10) ,AR,LAMDA,DELTA,NX,NY, 

* ALPHA, CLDSRD,ELLIP,XAC,CLA 
CHARACTERS IND 
IFdND.EQ. 'CONT' ) GOTO 100 
IF(IND.EQ. 'SELF' ) GOTO 110 

100 W = WHV1(WING(KL,10) ,WING(KL,12) ,WING(IJ,6) ,WING(IJ,8)) 

* - WHV1(WING(KL,10) ,WING(KL,12) ,WING(IJ,6) ,WING(IJ,9)) 

* - WHV1(WING(KL,10) ,WING(KL,12) ,WING(IJ,6) , -WING(IJ,8)) 

* + WHV1(WING(KL,10) ,WING(KL,12) ,WING(IJ,6) , -WING(IJ,9)) 

W = W*. 25/PI 

RETURN 

110 W = WHV2(WING(KL,11) ,WING(KL,12) ,WING(IJ,6) ,WING(IJ,8)) 

* - WHV2(WING(KL,11),WING(KL,12),WING(IJ,6),WING(IJ,9)) 

* - WHV2(WING(KL,11) ,WING(KL,12) ,WING(IJ,6) , -WING(IJ ,8)) 

* + WHV2(WING(KL,11),WING(KL,12),WING(IJ,6),-WING(IJ,9)) 

120 W = W*. 25/PI 

RETURN 

END 

-k:k:k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-!c-k-kirk-k-k-k-k-k-k-k-k-k->:-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k:k 
* * 

* FUNCTION WHV1(X1,Y1,X2,Y2) ^ 

ic ^ 

FUNCTION WHV1(X1,Y1,X2,Y2) 

IF (ABS(X1-X2) .LT. .0001) GOTO 100 

WHVl - (1.0+SQRT((X1-X2)**2 + (Y1-Y2)**2)/(X1-X2) )/(Yl-Y2) 
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RETURN 



100 WHVl = 1.0/(Y1-Y2) 

RETURN 

END 

'k k 

* FUNCTION WHV2(X1,Y1,X2,Y2) * 

k k 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

FUNCTION WHV2(X1,Y1,X2,Y2) 

IF (ABS(X1-X2) .LT. .0001) GOTO 100 

WHV2 = (1.0+(X1-X2)/SQRT((X1-X2)**2 + (Yl-Y2)**2) )/(Yl-Y2) 

RETURN 

100 WHV2 = 1.0/(Y1-Y2) 

RETURN 

END 



139 



MULTI 



•st -)5r -jS: * * -sV -st ^ -st -st * * * -st -3^ -jS: -A- -3^ "jS: * -3^ -st -jS: * * -st -5t •* * * -jS: ”5^ -st -jS: ■3^r *5^: 



•sSr * 

^ SUBROUTINE MULTI * 

* * 

* SUBROUTINE TO MULTIPLY THE SPECIFIED LOADING VECTOR, (1,17) * 

* BY THE COEFFICIENT MATRIX, AA( ), TO GET THE WING SHAPE * 

* VECTOR (1,19). * 

* * 

* VARIABLE NAMES ARE CONTAINED IN THE MAIN PROGRAM. * 

* ALL VARIABLES AND PARAMETERS ARE PASSED IN COMMON. * 

* THE ONLY VARIABLE CHANGED IS WING(I,19) * 

* * 



******************************* *********************** *********************** 
SUBROUTINE MULTI 
REAL LAMDA 
CHARACTER*3 ELLIP 

COMMON PI,B,WING(101,23) , SECTN(11 , 10) ,AR, LAMDA, DELTA, NX , NY, 

* ALPHA, CLDSRD, ELLIP, XAC,CLA 
COMMON/COF/A(101 , 101) , NEQNS 
COMMON/CAM/AA(101,101) , AS(101,101) 

DO 110 I = 1, NEQNS 
WING(I,19) = 0.0 
DO 100 J = 1, NEQNS 

WING(I,19) = WING(I,19)+AA(I,J)*WING(J.17) 

100 CONTINUE 



140 



110 CONTINUE 



RETURN 

END 
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MULT2 



*******:Ar**:A^******’A^')V'A'******>5t**:jt^**-5^:A-:A-***'3V****^**’3t>5^’5VV^**'A-**-jt**’3%'*’3V***'st’3t*^* 

•k k 

* SUBROUTINE MULT2 * 

* * 

* SUBROUTINE TO MULTIPLY THE FIAT PLATE LOAD VECTOR, (1,15) * 

* BY THE COEFFICIENT MATRIX, AS( ), TO GET THE DOWN WASH * 

* VECTOR (1,22). * 

* * 

* VARIABLE NAMES ARE CONTAINED IN THE MAIN PROGRAM. * 

* ALL VARIABLES AND PARAMETERS ARE PASSED IN COMMON. * 

* THE ONLY VARIABLE CHANGED IS WING(I,22) * 

* * 
•k-k-k-k-icirk-k-ic-k-k-kirJrk-k-k-k-k-k-k-k-k-k-kirk-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k'k'k 

SUBROUTINE MULT2 
REAL LAMDA 
GHARACTER*3 , ELLIP 

COMMON PI,B,WING(101,23) , SECTN(11 , 10) ,AR, LAMDA, DELTA, NX, NY, 

* ALPHA, CLDSRD,ELLIP,XAC,CLA 
COMMON/COF/A(101,101) , NEQNS 
COMMON/CAM/AA(101,101) , AS (101, 101) 

DO 110 I = 1, NEQNS 
WING(I,22) = 0.0 
DO 100 J = 1, NEQNS 

WING(I,22) = WING(I,22)+AS(I,J)*WING(J,15) 



100 



CONTINUE 



110 CONTINUE 



RETURN 

END 



143 



MULT3 



** *** '55: ** ■jSr* * Vr ************ -st ******** 'A’ "A” ')V * Vr -A- **** “A- * 



* * 

* SUBROUTINE MULT3 * 

* * 

* SUBROUTINE TO MULTIPLY THE SPECIFIED LOADING VECTOR, (1,17) * 

* BY THE COEFFICIENT MATRIX, AS( ), TO GET THE DOWN WASH * 

* VECTOR (1,23). * 

* * 

* VARIABLE NAMES ARE CONTAINED IN THE MAIN PROGRAM. * 

* ALL VARIABLES AND PARAMETERS ARE PASSED IN COMMON. * 

* THE ONLY VARIABLE CHANGED IS WING (I, 23) ^ 

* * 



*********************************************************************** 
SUBROUTINE MULT 3 
REAL LAMDA 
CHARACTER*3 ELLIP 

COMMON PI,B,WING(101,23) , SECTN(11 , 10) , AR, LAMDA, DELTA,NX ,NY, 

* ALPHA, CLDSRD, ELLIP, XAC,CLA 
COMMON/COF/A(101 , 101) , NEQNS 
COMMON/CAM/AA(101,101) , AS (101, 101) 

DO 110 I = 1, NEQNS 
WING(I,23) - 0.0 
DO 100 J = 1, NEQNS 

WING(I,23) = WING(I,23)+AS(I,J)*WING(J,17) 

100 CONTINUE 



144 



110 CONTINUE 



RETURN 

END 
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GRAP85 



icicicic'kic'k:kicic^'kic'kic'kic'kic'k'k'kic:k'kicicicicic'kicicicicicic'kicicicic^'kicic'k'k'k:k‘k:k'k'k'k‘k'k'k'k'k'k'k-k'k'k-kic'k'k-k'k 
•k k 

* PROGRAM GRAP85 * 

* * 

* THIS IS A PROGRAM WHICH GENERATES A MENU OF GRAPHIC OPTIONS * 

* AND WRITES FILES TO DISK WHICH ARE USED BY GRAPHER, A SOFTWARE * 

* PROGRAM FROM GOLDEN SOFTWARE, INC. BOX 281, GOLDEN, COLO., * 

* 80402. IF YOU DON'T HAVE A COPY OF THAT SOFTWARE, YOU MUST * 

* USE SOME OTHER METHOD OF PRESENTING THE DATA. * 

* * 

* DATA FILES USED FOR GRAPHIC OUTPUT ARE: * 

* 30 = X,Y DATA FILE * 

* 32 = GRAPH FILE * 

* * 

* * 

* ALL PARAMETERS ARE PASSED IN COMMON BLOCKS. * 

* VARIABLES USED ARE: * 

* * 

* ELLIP : 'YES' OR 'NO ', DEPENDING ON WHETHER THE SHAPE * 

* WAS COMPUTED * 

* TYPE : CHARACTER STRING 'FLAT ' OR 'CAMBERED' * 

* ITYPE : INTEGER, 1 IF FLAT, 2 IF CAMBERED * 

* TITLE : THE TITLE OF THE GRAPH, 6 POSSIBLE TITLES ARE USED * 

* ISELEC : INTEGER RESPONSE FROM THE GRAPHIC SCREEN * 

* XMIN : MIN VALUE ON HORIZONTAL AXIS FOR PLOT OF A SECTION * 
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■k 


XMAX 


k 


YMIN 


k 


YMAX 


k 


VMAX 


k 


VMIN 


k 


VERT 


k 


lOFFS 


k 




k 


ISECT 


k 


ESC 


k 


HSTART 


k 




k 


VINC 


k 


XINC 


k 


YINC 


k 


PLE 


k 


PTE 


k 


XLE 


k 


XTE 


k 


AXIS 


k 




k 


TITLE 


k 




k 


TITLE2 


k 




k 


TYPE 



MAX VALUE ON HORIZONTAL AXIS FOR PLOT OF A SECTION * 
MIN VALUE ON HORIZONTAL AXIS FOR PLOT OF A SPAN * 
MAX VALUE ON HORIZONTAL AXIS FOR PLOT OF A SPAN * 
MAX VALUE ON VERTICAL AXIS OF PLOT * 
MIN VALUE ON VERTICAL AXIS OF PLOT * 
A TEMPORARY VARIABLE USED TO HOLD A VALUE. * 
AN OFFSET USED TO LOCATE THE CORRECT COLUMN IN THE * 
WING MATRIX * 
SPANWISE SECTION TO BE PLOTTED * 
THE ESCAPE CHARACTER (27) * 
POSITION ON THE GRAPH WHERE THE HORIZONAL AXIS * 
STARTS * 
VERTICLE AXIS INCREMENT IN UNITS * 
HORIZONTAL AXIS INCREMENT WHEN PLOTTING CHORD * 
HORIZONTAL AXIS INCREMENT WHEN PLOTTING SPAN * 
PLOT LEADING EDGE POSITION IN INCHES * 
PLOT TRAILING EDGE POSITION IN INCHES * 
LEADING EDGE POSITION IN PLANFORM UNITS * 
TRAILING EDGE POSITION IN PLANFORM UNITS * 
CHARACTER ARRAY WHICH HOLDS THE TITLES FOR THE * 
HORIZONTAL AXIS. * 
CHARACTER ARRAY WHICH HOLDS THE TITLES FOR THE * 
VERTICAL AXIS. * 
CHARACTER ARRAY WHICH HOLDS PART OF THE GRAPH * 
TITLE . * 
CHARACTER ARRAY WHICH HOLDS PART OF THE GRAPH * 
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* TITLE * 

* COEF : CHARACTER ARRAY WHICH HOLDS PART OF THE GRAPH * 

* TITLE. * 

* COEFF : REAL ARRAY WHICH HOLDS THE LIFT, DRAG OR MOMENT * 

* ASSOCIATED WITH A PARTICULAR GRAPH. * 

* * 

* SUBROUTINES CALLED ARE: * 

* CDAT85 , FORTRAN, GENERAL DATA ENTRY * 

^ DELAY , FORTRAN, SLOW DOWN THE PROGRAM SO ERROR MESSAGES * 

^ CAN BE SEEN BEFORE BEING ERASED, ^ 

* SDIG , FORTRAN, USED TO FIND THE NUMBER OF SIGNIFICANT ^ 

^ DIGITS IN THE AXIS LIMITS AND MAKE AXIS LIMITS * 

* MORE UNIFORM. * 

* RNDUP , FORTRAN, USED TO ROUND THE AXIS LIMITS TO THE NEXT * 

* WHOLE DIGIT. * 

■A- -k 

k k 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

PROGRAM GRAP85 
REAL LAMDA 
CHARACTER* 3 ELLIP 
CHARACTER*8 TYPE(2) 

CHARACTER*20 TITLE(6) 

CHARACTER*8 TITLE2(2,2) 



CHARACTER*! ESC 



CHARACTER*23 AXIS (2) 



CHARACTER*? COEF(2,3) 

REAL C0EFF(2,4) 

COMMON PI,B,WING(101,23) , SECTN(11 , 10) , AR, LAMDA,DELTA,NX,NY, 

* ALPHA, CLDSRD,ELLIP,XAC,CLA 
COMMON/COF/A(101,101) , NEQNS 
COMMON/CAM/AA(101,101) , AS (101, 101) 

COMMON/FILS/I IN , lOUT , JOUT , KOUT , LOUT , MOUT , NOUT 

* "k 

* DEFINE THE CHARACTER STRINGS NEEDED FOR GRAPH TITLES * 

* * 

TYPE(1)=' FLAT' 

TYPE (2)=' CAMBERED' 

TITLE(1)= 'd(CL)/dy 

TITLE(2)= 'd(CD)/dy 

TITLE(3)= 'd(CMAC)/dy 

TITLE(4)= 'TWIST IN DEGREES 

TITLE(5)= 'DELTA GP 

TITLE(6)= 'WING HEIGHT COORD Z ' 

TITLE2(1,1)= ' , AOA = ' 

TITLE2(1,2)= ' deg 
TITLE2(2,1)= ' , CLi = ' 

TITLE2(2,2)= ' 

AXIS(1)= 'SPAN-WISE COORDINATE Y ' 
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AXIS (2)- 'CHORD -WISE COORDINATE X' 



COEF(l,l)='CL = ' 

COEF(l,2)='CD = ' 

COEF(l,3)='CMAC = ' 

COEF(2,l)='CLi = ' 

COEF(2,2)='CD = ' 

COEF(2,3)='CMAC = ' 

* * 

* ASSIGN DEFAULT VALUES TO THE VARIABLES * 

* * 



*************************************************************************************** 



IMOD = 


1 


AR = 


o 

CM 


LAKDA = 


1.0 


ADELTA= 


1.0 


NX = 


4 


NY « 


4 


ALPHA = 


5.0 


ELLIP = 


'YES' 


CLDSRD= 


1.0 


ITYPE - 


1 


ISECT = 


1 


XAC = 


0.0 


ESC = CHAR(27) 


PI = 2. 


0*ACOS(0 
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***************************************■***■*********'******************'** 



* SET INPUT/OUTPUT VARIABLES AND OPEN THE NECESSARY FILES 



* 

* 



IIN - 5 
lOUT = 6 

JOUT - 10 

KOUT =11 
LOUT =12 
LDATl = 13 
LDAT2 = 14 
OPEN(UNIT = IIN) 

OPEN(UNIT = lOUT) 

OPEN(UNIT = JOUT, FILE= 'WING. MAT ' ) 

OPEN(UNIT = KOUT, FILE=' SECTION. MAT' ) 

OPEN(UNIT = LOUT, FILE=' PARAMS .DAT' ) 

OPEN (UNIT = LDATl, FILE=' FLAT.DAT' ) 

OPEN(UNIT = LDAT2, FILE=' CAMBER.DAT' ) 

10 FORMAT (I6/3(F6.2,/) ,2(13/) ,F6.2/A/F6.2/I2/I2/F6.2) 

20 FORMAT (/2 (lOX, F12 . 6/)/10X, F12 . 6) 

'A' -sSr -st ^ ')Sr 'sSr -:Sr *5^ *3t * -A- -A- -3^ *55: 'A' 'A' -5^ '*■ ^ 



* FORMAT STATEMENTS ARE FOR THE INPUT TO THE GRAPHER PROGRAM. * 

* OTHER GRAPHICS PROGRAMS WILL REQUIRE DIFFERENT FORMAT STMTS. * 

* * 
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^ •A' ^ -it •sV ^ ^ *35: ’A’ ^ 'A’ * -it -iJr ^ "jV "jS: ^ ^ •A' ^ -iSr "sSr ^ ^ "sSr * ^ ^ * 



32 FORMAT ( 

* '1236'/ 

* '1 2 1 0 O'/ 

* 'PI'/ 

* '65 66 78 "NO "O'/ 

* '"NO" "SOLID" 1.500e-001 1'/ 

* '"YES" 41 l.OOOe-OOl 1 1'/ 

* '78 9.900e+028 9.900e+028 O.OOOe+000 "DEFAULT" l.OOOe-OOl 1'/ 

* '"SOLID" 5 1.500e-001 O.OOOe+000 9.900e+029 100 2.000e+000 1') 

33 FORMAT ( 

* '1236'/ 

* '1 2 3 0 1'/ 

* 'PI'/ 

* '65 66 78 "NO " 0'/ 

* '"NO" "SOLID" 1.500e-001 1'/ 

* '"YES" 41 l.OOOe-OOl 1 1'/ 

* '78 9.900e+028 9.900e+028 O.OOOe+000 "DEFAULT" l.OOOe-OOl 1'/ 

* '"SOLID" 5 1.500e-001 O.OOOe+000 9.900e+029 100 2.000e+000 1') 

34 FORMAT ( 

* 'HORIZ'/ 

* '1.750e+000 '.E10.3,' 6.000e+000 88'/ 

* E10.3,' '.E10.3,' 9.900e+028 11'/ 

* 'O.OOOe+000 '.E10.3,' 1.500e-001 1 1'/ 

* '1 1 1 '/ 

* '1 2.750e+000 4.000e-001 9.900e+028 1.800e-001'/ 
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* ' "DEFAULT" "DEFAULT" " ' , A, ' " ' ) 



36 FORMAT ( 

* 'Y-AXIS'/ 

* '1.750e+000 l.OOOe+000 5.000e+000 89'/ 

* E10.3,' '.E10.3,' 9.900e+028 1 1'/ 

* '2.700e+002 '.E10.3,' 1.500e-001 1 1'/ 

* '1 ' , 11 , ' 1 '/ 

* '1 9.900e+028 O.OOOe+000 9.900e+028 1.800e-001'/ 

* '"DEFAULT" "DEFAULT" "',A,'"') 

37 FORMAT ( 

* 'PI'/ 

* '"DEFAULT" 1'/ 

* '2.000e+000 7.000e+000 O.OOOe+000 1.700e-001'/ 

* ' 2 '/ 

* '0 13 "' ,A, ' WING"'/ 

* '1 12 "CLi = ' ,F5.2, '"') 

38 FORMAT ( 

* 'PI'/ 

* '"DEFAULT" 1'/ 

* '2.000e+000 7.000e+000 O.OOOe+OOO 1.700e-001'/ 

* ' 2 '/ 

* '0 13 "' ,A, ' WING"'/ 

* '1 34 "' ,A,F6.4,A,F5.2,A'"') 

39 FORMAT ( 

* 'PI'/ 

* '"DEFAULT" 1'/ 
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* '2.000e+000 7.000e+000 O.OOOe+000 1.700e-001'/ 

* ' 2 '/ 

* '0 26 WING, SECTION ',13,'"'/ 

* '1 26 "vertical scale exaggerated"') 

40 FORMAT ( 

* 'PI'/ 

* '"DEFAULT" 1'/ 

* '2.000e+000 7.000e+000 O.OOOe+000 1.700e-001'/ 

* ' 1 '/ 

* '0 26 "',A,' WING, SECTION ',13,'"') 

41 FORMAT ( 

* 'LE'/ 

* '"DEFAULT" 1'/ 

* E10.3,' 6.000e+000 O.OOOe+OOO l.OOOe-OOl'/ 

* ' 1 '/ 

* '0 4 "L.E."') 

42 FORMAT ( 

* 'TE'/ 

* '"DEFAULT" 1'/ 

* E10.3,' 6.000e+000 O.OOOe+OOO l.OOOe-001'/ 

* ' 1 '/ 

* '0 4 "T.E."') 

43 FORMAT ( 

* 'PI'/ 

* ' 2 '/ 

* E10.3,' l.OOOe+000 ',E10.3,' 6.000e+000 1 2 l.OOOe-001'/ 



1 “=' 



50 



* E10.3,' l.OOOe+000 '.ElO.3,' 6.000e+000 1 2 l.OOOe-001') 
FORMAT (2(F10.6,1X)) 



60 FORMAT (23F10.6) 

70 FORMAT (10F10.6) 

St*******************************************'********'***'*'*******'*'****'**'* 
* * 

* THESE FORMAT STATEMENTS ARE FOR THE FLASHUP WINDOWS MENUS * 

* * 
***************************'!lr*******-****-********')SriV*****-***********'****** 

82 FORMAT (' ' ,A,A,I3) 

84 FORMAT (' ' ,A,A,A) 

'kic-k-k'k-k-kic‘k’k'k":k-k’k'k-k‘k'k'k'k'k-k-k’kii'k'k‘k‘k’k'k‘k'k-k'k'kic-k'k’k'k'k‘k 

■k -k 

* READ THE LAST SET OF PARAMETERS, IF UNCHANGED, READ WING AND * 

* SECTN DATA , IMOD IS THE FLAG TO TELL IF THE PARAMETERS HAVE * 

* BEEN CHANGED. 1 = CHANGE, 0 = NO CHANGE. * 

* * 
*********************************************************************** 

READ (LOUT, 10) IMOD , AR , LAMDA , ADELTA , NX , NY , ALPHA , 

* ELLIP,CLDSRD,ITYPE,ISECT,XAC 
DELTA = ADELTA*PI/180. 

NEQNS = NX*NY 
B = AR*(1.0+LAMDA)/2.0 
IF (IMOD.EQ.O) THEN 

READ (JOUT,60) ( (WING(I ,J ) ,J=1 , 23) , 1=1 , NEQNS) 

READ (KOUT,70) ( (SECTN(I , J ) , J=1 , 10) , 1=1 ,NY) 
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READ (LDAT1.20) (C0EFF(1 , J) , J=1 . 3) 

READ (LDAT2.20) (C0EFF(2 , J) , J=1 , 3) 

ELSE 

WRITE (lOUT,*) ' CHANGES WERE MADE TO THE PLANFORM. GO BACK' 
WRITE (lOUT,*) ' AND RUN THE LOAD/SHAPE PROGRAM AGAIN.' 

CALL DELAY (6) 

GO TO 1000 
ENDIF 

C0EFF(1,4) = ALPHA 
COEFF(2,4) = CLDSRD 
CLOSE (JOUT) 

CLOSE (KOUT) 

CLOSE (LDATl) 

CLOSE (LDAT2) 

*********************************************************************** 
* * 

* SET MINIMUM VALUE OF X. THIS X IS THE CHORDWISE COORDINATE * 

* * 
*********************************************************************** 

XMIN =0.0 

•k k 

* COMPUTE THE CORNER OF THE WING, TO GET MAX VALUES FOR THE AXIS * 

* THIS IS NEEDED IN CASE THE WING IS SWEPT OR HIGHLY TAPERED. * 

* * 
************************************************************************ 



156 



XT = ( B*TAN( DELTA) /2.0) + LAMDA 



•k ‘k 

* NOW CONVERT SPAN AND CHORD AXIS LIMITS INTO "ROUNDED” UNITS * 

* AT THE SAME TIME, COMPUTE THE INCREMENTS IN TIC MARKS FOR THE * 

* GRAPHS . * 

k k 



kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 



IF (XT.GT.l.) THEN 

CALL RNDUP(XT.XMAX) 

ELSE 

XMAX = 1.0 
ENDIF 

YMIN =0.0 

CALL RNDUP(B/2.0,YMAX) 

XING = (XMAX-XMIN)/5.0 
YINC = (YMAX-YMIN)/5.0 

************ *******************'Jr**i(r***-****'it***')V')t**'>lr**V<r'jt'j!r**'ilr**'5lf'!t*'>t')<:**'jlr*'>!r 
* * 

* CLEAR THE SCREEN AND WRITE THE GRAPHIC MENU. * 

* * 
*********************************************************************** 
110 WRITE (lOUT,*) ESC,'[2J' 

WRITE (lOUT,*) ESC, ' [00;00H' 

WRITE (lOUT,*) '~C=ALL/' 
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WRITE (TOUT,*) ' W=GRAPH/' 

WRITE (IOUT.84) ESC , ' [08 ; 69H' , TYPE(ITYPE) 
WRITE (IOUT.82) ESC, ' [09 ; 69H',ISECT 



'k k 

* READ THE OPTION AND SEND TO THE APPROPRIATE AREA. * 

* OPTIONS ARE: * 

* 1, SETS THE TYPE OF WING, FLAT OR ELLIPTIC. * 

* 2, SETS THE WING SECTION TO BE PLOTTED * 

* 3, RETURNS TO THE MAIN MENU * 

* 4-7 PLOT THE LIFT, DRAG, MOMENT OR TWIST VS SPAN * 

* 8 PLOT DELTA CP VERSUS CHORD * 

* 9 PLOTS THE WING SHAPE VERSUS CHORD, WHICH GIVES * 

* AN ELLIPTIC LIFT DISTRIBUTION. * 

* * 



READ (IIN,*) ISELEC 
WRITE (lOUT,*) ESC,'[2J' 

WRITE (lOUT,*) ESC, ' [00;00H' 

************************************************************************************** 
* * 

* THIS IS OPTION 1, SET TYPE OF WING. ^ 

* * 
kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk'kkkk'kkkkkkk 

IF (ISELEC. EQ.l) THEN 

IF (ELLIP.EQ. 'NO ') THEN 



158 



WRITE (lOUT,*) ' ELLIPTIC DATA WAS NOT COMPUTED' 

CALL DELAY(6) 

GO TO no 

ELSE 

ENDIF 

IF (ITYPE.EQ.l) THEN 
ITYPE = 2 

ELSE 

ITYPE = 1 

ENDIF 

GO TO no 

ELSE 

ENDIF 

*********************************************************************** 
* * 

* THIS IS OPTION 2, SET SECTION NUMBER * 

•k k 



kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

IF (ISELEC.EQ.2) THEN 

CALL CDIT85(1,NY,ISECT, 

* 'WING SECTION NUMBER ',IMOD) 

IMOD = 0 

GO TO no 

ELSE 

ENDIF 



kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 
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'k -k 

^ THIS IS OPTION 3, RETURN TO MAIN PROGRAM. BEFORE LEAVING, WRITE * 

* THE FLASHUP WINDOW FOR THE MAIN MENU. * 

* 'k 



ic'kk'ki^'kk'kk'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'kk'k'k'kk'k'k'k'k'k'k-k'k'k'kk'k-k'k-k'k'k'k'k'k'k'k'k'k-k'k'k 

IF (ISELEC.EQ.3) THEN 
GO TO 1000 
ELSE 
END IF 

IF ((ISELEC.LT.l) .OR. (ISELEC.GT.9)) THEN 

WRITE (lOUT,*) ' THAT IS NOT A VALID SELEGTION' 

CALL DELAY(6) 

GO TO 110 
ELSE 
END IF 

************************************************************************************ 



* * 

* THESE ARE THE LIFT, DRAG , MOMENT AND TWIST VS SPAN PLOTS. * 

'k -k 

'k'kic'kic'k'k'kk'kk'k'kk'k'k'k'k'k'k'k'k'k'k'kic'k'k'k'kk'k'k'k'kicic'k'k'kkkiic'ki^ic'k'k'k'k'k'kk'k'k'k'k'k'k'k'k'kk'^ 

IF (ISELEC.GT.7) GO TO 310 

'k’k'k'k'k'k'k'k-k'k-k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k:k'k'k'k'k'k'k'k'k'k'k'k'kk'k'k'k 
'k :k 

* OFFSET IS ESTABLISHED, IT IS AN EASY WAY TO USE THE ISELEC * 

* TO FIND THE CORRECT COLUMN OF THE SECTN MATRIX. * 

* SINCE THERE IS NO TWIST ON A FLAT WING, SET THE WING TYPE TO * 






n 






CAMBERED IF THE TWIST PLOT IS SELECTED. 



* 

* 






IF (ISELEC.EQ.7) ITYPE = 2 
IF (ITYPE. EQ.l) THEN 
lOFFS = -2 
ELSE 

lOFFS = 3 
ENDIF 

icic'k-kic-icic-kic-icic-kic-ic-ic-kic-ic-ic-kicic-ic-k-^-ic-k-k-^-icic-k-^-icic-k-k-ic-irk-kic-k-k-^icic-k-^-ic-ic-k-^-irk-k^-k-k'irkic'^icJric^lrA-kie-k 
* * 



* SET UP THE DATA FILE FOR THE GRAPHIC PLOT. * 

* ■* 



*********************■*■*■*■■*■■**■*■*■*■*■*■*■*■!(:■*■*■*****■*■*■*■*■***■!<:■)!:**■!(:■*• *'sV'A:****'*’'****'sV>!r** 

J = 1 

VMAX = -100.0 
VMIN = 100.0 

OPEN (UNIT = 30, FILE ='P1.DAT') 

IJ = J 

DO 300 I = 1,NY 

IJ = (I-1)*NX + J 

VERT = SECTN(I,IOFFS+ISELEC) 

IF ( VMAX. LT. VERT) VMAX = VERT 
IF ( VMIN. GT. VERT) VMIN = VERT 
WRITE (30,50) WING(IJ,5), VERT 
300 CONTINUE 
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CLOSE (30) 






'k k 

* SET UP THE Pl.GRF FILE FOR GRAPHER. IT HAS ALL THE GRAPH * 

* FORMATS IN IT. * 

* * 



***************** ****************************i<r*****************'*****'*''*'* 

IF (VMIN.GT.O.O) VMIN =0.0 
IF (VMAX.LT.0.0) VMAX = 0.0 
CALL RNDUP (VMAX, VMAX) 

CALL RNDUP (VMIN, VMIN) 

CALL SDIG( VMIN, VMIN, VMAX, VMAX, IDIG) 

VINC = (VMAX- VMIN) /5.0 

HSTART =(VMAX-6*VMIN)/(VMAX-VMIN) 

OPEN(UNIT = 32, FILE ='P1.GRF') 

WRITE (32,32) 

WRITE (32,34) HSTART, YMIN,YMAX,YINC, AXIS (1) 

WRITE (32,36) VMIN,VMAX,VINC, IDIG,TITLE(ISELEC-3) 

IF (ISELEC.EQ.7) THEN 

WRITE (32,37) TYPE(ITYPE) ,COEFF(2 ,4) 

ELSE 

WRITE (32,38) TYPE(ITYPE) , COEF(ITYPE , ISELEC-3) , 

* COEFF(ITYPE, ISELEC-3) , TITLE2(ITYPE, 1) , C0EFF(ITYPE,4) , 

* TITLE2(ITYPE,2) 

ENDIF 

CLOSE (32) 
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GO TO 1000 



■:t-3V'3Sr'A'Vt':V-A':Ar'it'3V'3Sr:A'-3VVr:5Sr:A-'3SrVr'3Sr‘3V'5Sr'3Sr'35:-3t:;Sr:A':5V:>V-:Sr'5Sr‘A''5V^ 

* * 

^ THESE ARE THE DELTA CP PLOTS * 

* * 

•:k'k-kic'kic-k'kic':kic-:k'k'k-k':k'k'kicic'k'k-:k-:k'k'k'kicic-:k-k'k'k'kicicic'kic'kicic'kici(icic'kicicicicicicicic-ki^ 



310 CONTINUE 

IF (ISELEC.EQ.9) GO TO 400 



Vr ic 

* OFFSET IS SET AGAIN. 

* 

'k:k'k:k'k:k'k'k'kic'k-k':kic':k-kic-:kic-:kicicic‘:kic-:kicicicicicicic-:kic-:kii^icicic^i(icicicicic:ki^iiifi(icicicicis-:ki^ic‘:kicicicicicic-:k-:kic 



IF (ITYPE.EQ.l) THEN 
lOFFS = 8 
ELSE 

lOFFS - 10 
ENDIF 

********************************************************************************* 
* * 



* SET UP THE PI. DAT FILE FOR GRAPHER 

* 



* 

* 



************************************■*•***■*•■*•**■*•■*•*■*■*•***■*•***•*■**■*•**■*■*•**■*•**•*•■*• 
I = ISECT 
VMAX = -100.0 
VMIN = 100.0 
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OPEN (UNIT = 30, FILE ='P1.DAT') 
IJ = (I-1)*NX + 1 



'k 

ic 

ic 

* 

■k 



XLE IS THE VALUE OF THE LEADING EDGE. SINCE THE 
FOR THE CAMBERED WING IS ZERO AT THE LEADING AND 
EDGE, IT IS POSSIBLE TO EXTEND THE DATA TO THOSE 





•A 


DELTA CP 


■A 


TRAILING 


A- 


POINTS . 


A- 




A 



^k-k'kic^kic'kicici^'k'k’k'k-k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k-k'k-k'k’k'k'k'k'k'k'k'k'k-k'k'k'k'k-k'k'k'k-k'k'k'k 

XLE = WING(IJ,4)-WING(IJ,13)/2. 

IF(ITYPE.EQ.2) THEN 

WRITE (30,50) XLE, 0.0 
ELSE 
END IF 

'k:k:k-k'kic:k-k-k:k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k-k'k'k'k'k'k'k'k'k'k'k'k'k:k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 



^k 



'k 



* GET THE MAX AND MIN VALUES OF THE PARAMETER BEING PLOTTED ^ 

* SO ALL GRAPHS WITH THAT PARAMETER WILL HAVE THE SAME VERTICAL 

^ SCALE. ^ 



'k 



'k 



* •3^ ^ ^ 'A- -55: 'A- "sS: -sSr -jS: •:*: ”3!^ * -A- -55: •jSr -jS: -55: -55: -)5r -A- *3^ -sJ: ”3^ "st -55: 

DO 315 IJ = l.NEQNS 

VERT = WING(IJ,IOFFS+ISELEC) 
IF (VMAX.LT.VERT) VMAX = VERT 
IF (VMIN.GT.VERT) VMIN = VERT 
315 CONTINUE 



16A 



•>\r ->V 'A' "sV ■>5r - jV -jV *5^ •>5r - jV • sV -jV -st - jV - jt -jSr 



^ SET UP THE PI. DAT FILE FOR GRAPHER ^ 

■k 'k 

DO 320 J = 1,NX 

IJ = (I-1)*NX +J 

VERT = WING(IJ ,IOFFS+ISELEC) 

WRITE (30,50) WING(IJ,11), VERT 
320 CONTINUE 
IJ = I*NX 

■:*: -st *:Sr "st •>V -st •>V "sV •>V -jV "sV -sV •>V ■>V ->V '5^ •>V’ ^ ->V •j^r "jV •>V’ •>V -^Sr ->V "st •>V •>V •sV ■5^ •>V’ ^ -sV ->5r 



* XLE IS COMPUTED FOR THE REASONS STATED ABOVE. ^ 

^ 'k 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk'k'kkkkkkkkkkkkkkk'kk'kkkkkk'kk'kkkkk 



XTE = WING(IJ,4)+WING(IJ,13)/2. 

IF(ITYPE.EQ.2) THEN 

WRITE (30,50) XTE, 0.0 
ELSE 
ENDIF 
CLOSE (30) 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 
k k 



k 

k 



COMPUTE AND STORE THE VALUES FOR THE Pl.GRF FILE USED BY GRAPHER* 

* 
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'k'k'kici<:ic'k‘^rk:k:k:k:k:k-:k:k:k'k'k'^(:k'k:k:k:k'k'^'k'kic'k'kicic'kic'k-k'k'k'kic'k'k'kicic'ki^'k'kic'kic'k'kic'k'k'k 



IF (VMIN.GT.0.0) VMIN =0.0 
IF (VMAX.LT.0.0) VMAX - 0.0 
CALL RNDUP( VMAX, VMAX) 

CALL RNDUP( VMIN, VMIN) 

CALL SDIG (VMIN, VMIN, VMAX, VMAX, IDIG) 

VINC = (VMAX- VMIN) /5.0 

HSTART = (VMAX- 6*VMIN)/ (VMAX- VMIN) 

PLE = (6.*XLE/(XMAX-XMIN))+1.75 
PTE = (6.*XTE/(XMAX-XMIN))+1.75 
0PEN(UNIT = 32, FILE ='Pl.GRF') 

WRITE (32,33) 

WRITE (32,34) HSTART, XMIN,XMAX,XINC,AXIS(2) 

WRITE (32,36) VMIN, VMAX, VINC, IDIG, TITLE(ISELEC-3) 

WRITE (32,40) TYPE(ITYPE) , ISECT 
WRITE (32,41) PLE 
WRITE (32,42) PTE 
WRITE (32,43) PLE , PLE, PTE, PTE 
CLOSE (32) 

GO TO 1000 

*********************************************************************** 



* THIS IS THE SHAPE PLOT. 



* 

* 



*********************************************************************** 



400 IF (ELLIP.EQ. 'NO ') THEN 









WRITE (lOUT,*) ' ELLIPTIC SHAPE DATA WAS NOT COMPUTED' 

CALL DELAY(6) 

GO TO 110 

ELSE 

ENDIF 

**************************************'*******************'****★**■*•****** 
* * 

* SET UP THE PI. DAT FILE FOR THE WING SHAPE. WE KNOW THAT THE * 

* LEADING EDGE STARTS AT ZERO HEIGHT, SO THAT STARTING VALUE IS * 

* ADDED TO THE DATA, IN ORDER TO GET THE MAXIMUM INFORMATION * 

* INCLUDED IN THE GRAPH * 

* XLE AND XTE ARE ALSO NEEDED SO THAT VERTICAL DOTTED LINES * 

* SHOWING THEIR POSITION CAN BE ADDED TO THE GRAPH. * 

* * 
********************************************************* -A-****************** 

I = ISECT 
VMAX = -100.0 
VMIN = 100.0 

OPEN (UNIT = 30, FILE ='P1.DAT') 

IJ = (I-1)*NX + 1 

XLE = WING(IJ,4)-WING(IJ,13)/2. 

WRITE (30,50) XLE, 0.0 
IJ - I*NX 

XTE - WING(IJ,4)+WING(IJ,13)/2. 

DO 405 IJ = 1,NEQNS 
VERT = WING ( I J, 20) 
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IF (VMAX.LT.VERT) VMAX = VERT 
IF (VMIN.GT.VERT) VMIN = VERT 
405 CONTINUE 

DO 410 J = 1,NX 
IJ = (I-1)*NX+J 
VERT = WING (IJ, 20) 

WRITE (30,50) WING(IJ,4)+WING(IJ,13)/2.0, VERT 
410 CONTINUE 

CLOSE (30) 

IF (VMIN. GT. 0.0) VMIN =0.0 
IF (VMAX. LT. 0.0) VMAX =0.0 
CALL RNDUP (VMAX, VMAX) 

CALL RNDUP (VMIN, VMIN) 

CALL SDIG (VMIN, VMIN, VMAX, VMAX, IDIG) 

VINC = (VMAX- VMIN) /5.0 

HSTART = (VMAX- 6*VMIN)/ (VMAX- VMIN) 

PLE = (6.*XLE/(XMAX-XMIN))+1.75 
PTE = (6.*XTE/(XMAX-XMIN))+1.75 
OPEN(UNIT = 32, FILE ='P1.GRF') 

WRITE (32,33) 

WRITE (32,34) HSTART,XMIN,XMAX,XINC,AXIS(2) 

WRITE (32,36) VMIN , VMAX, VINC , IDIG ,TITLE(ISELEC- 3) 

WRITE (32,39) TYPE(2) , ISECT 

WRITE (32,41) PLE 

WRITE (32,42) PTE 

WRITE (32,43) PLE , PLE , PTE , PTE 
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CLOSE (32) 

1000 REWIND (LOUT) 

***************************************************************************** 
* * 

* RE-WRITE THE PARAMETER DATA FILE, SINCE IT KEEPS TRACK OF THE * 

* LAST WING TYPE AND SECTION PLOTTED. * 

* * 
*******************************************************■***'*******■****** 

WRITE (LOUT, 10) IMOD , AR , LAMDA , ADELTA , NX , NY , ALPHA , 

*ELLIP , CLDSRD , ITYPE , I SECT , XAC 
CLOSE (LOUT) 

ic * 

^ BEFORE LEAVING THE PROGRAM, RE-WRITE THE MAIN PROGRAM MENU. * 

^ * 

WRITE (lOUT,*) ESC,'[2J' 

WRITE (lOUT,*) ESC, ' [00;00H' 

WRITE (lOUT,*) '‘C=ALL/' 

WRITE (lOUT,*) ''W=MAIN/' 

END 

$ INCLUDE: 'CDIT85.FOR' 

$ INCLUDE : ' DELAY . FOR ' 

$INCLUDE: 'RNDUP.FOR' 

$INCLUDE: 'SDIG.FOR' 
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DELAY 



•:k'k'k'k'k'kic-kicic'k^'ki^^:k'krk'k^'k-k'k'k'k'k:k^:k:k'k'k:k’kic'k'k-k'k:k'k'k'k'k'k'k'k'k'k'k'k:k'k'k'k'k'k:k'k-k'k:k'k:k'k'k'k'k'k'k'k 



* 










k 


* 


SUBROUTINE 


, DELAY(INTERV) 






k 


'k 










k 


* 


CALLED TO 


SLOW THE SCREEN 


DOWN AFTER PRINTING 


AN ERROR MESSAGE. 


k 


'k 


VARIABLES 


USED ARE: 






k 


-k 


INTERV 


: NUMBER OF SECONDS FOR DELAY MUST 


BE LESS THAN 59, 


k 


•k 




INPUT. 






k 


k 


IHR 


: TIME OF DAY 


IN HOURS 




k 


k 


IMIN 


: TIME OF DAY 


IN MIN 




k 


k 


ISECl 


: TIME OF DAY 


IN SECONDS 




k 


k 


ISEC2 


: TIME OF DAY 


IN SECONDS 




k 


k 


IHUND 


: TIME OF DAY 


IN HUNDREDTHS OF SECONDS. 


k 


k 










k 



'kic'k-k'k'k'k-k^ic'k:k-kic'kicic'kic'k'kicic':kic'k'k’k’k:k'kic'k'k:k'k’k'k'k'k'kic'k-k'k’k'k9c'k'k'k^ic'k':k'k'k'k'k’k-k'k'k:k:ki^^ 

SUBROUTINE DELAY(INTERV) 

INTEGER*2 IHR, IMIN, ISECl, ISEC2, IHUND, INTER 
INTER = INTERV 

CALL GETTIMC IHR, IMIN, ISECl, IHUND) 

100 CALL GETTIM(IHR, IMIN, ISEC2, IHUND) 

IF ( (ISEC2- ISECl ) .LT. INTER) GOTO 100 

RETURN 

END 
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CDIT85 



ici(ic'^'kicic'ki<'k-ki^i<:-^icicic^i(icicic'k'kic'kicic'ki:i:'kic'k'k:k'ki^i(ic‘k'ki('k'k'k'k'k’k'k'k'k'k-k'k'k'k'k'k'k':k'k'k'k'k'k'k'k'k'k'k 



'k 










k 


SUBROUTINE 


CDIT8 5 ( LOW , HI GH , VALUE , NAME , IMOD ) 




k 


k 


INTEGER VERSION OF CDAT85 




k 


k 


HANDLE INPUT OF PARAMETER DATA FOR THE PROGRAM. INCLUDES SOME 


k 


k 


ERROR CHECKING. 




k 


k 








k 


k 


VARIABLES 


USED ARE: 




k 


k 


LOW 


: MINIMUM VALUE THAT WILL BE ACCEPTED, 


INPUT 


k 


k 


HIGH 


: MAXIMUM VALUE THAT WILL BE ACCEPTED, 


INPUT 


k 


k 


VALUE 


: VALUE READ FROM SCREEN, OUTPUT 




k 


k 


NAME 


: SCREEN PROMPT STRING, INPUT 




k 


k 


IMOD 


: FLAG USED IN ANOTHER PROGRAM, SET TO 


1 ANYTIME 


k 


k 




CDAT85 IS CALLED., OUTPUT 







k 


ESC 


: ESCAPE CHARACTER, (27) 


k 


k 






k 



'k'k'kic'kic'k'k'k:ki^i('k'^'kic'k'k'k'k'k'k'k'k^'k'k'kic'ki('ki(icicicicici(ic'kif:'kici(i^ic'ki(i(ic'kicic'ki(i(i<i(ici(ici(‘‘k'kic‘k'k'k'kic 

SUBROUTINE CDIT85 (LOW , HIGH , VALUE , NAME , IMOD) 

COMMON/FILS /I IN , I OUT , JOUT , KOUT , LOUT , MOUT , NOUT 

INTEGER LOW, HIGH, VALUE 

CHARACTER*35 NAME 

CHARACTER*! ESC 

ESC = CHAR(27) 

10 FORMATC' ENTER THE ' ,A35) 
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20 



FORMAT(' THE MINIMUM VALUE FOR ',A35,/' IS ',13) 



30 FORMATC THE MAXIMUM VALUE FOR '.A35,/' IS ',13) 

^ 'A' ‘:iV -jSr -sSr ^ •:Sr •:Sr ■A* •:Ar ■jV •3^r -jV - sSr •:Sr •A' •A' 



■A A 

* CLEAR SCREEN, POSITION CURSOR, AND * 

* SET IMOD FLAG AND WRITE PROMPT TO SCREEN. ^ 

A A 



AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

100 WRITE (lOUT,*) ESC,'[2J' 

WRITE (lOUT,*) ESC,'[10;0H' 

IMOD = 1 

WRITE (lOUT, 10) NAME 

AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 



A A 

* READ SCREEN VALUE. ON ERROR OR END OF DATA SET, GO BACK TO * 

* PROMPT . ^ 

A A 



AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

READ(IIN,*,ERR = 100, END = 100) VALUE 

AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 



A A 

* IF VALUES ARE OUT OF RANGE , PRINT ERROR MSG AND RETURN TO * 

* PROMPT . ^ 

A A 



AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 



1 ^- 






IF (VALUE.lt. LOW) THEN 



WRITE (IOUT.20) NAME, LOW 
CALL DELAY (6) 

GO TO 100 
ELSE 

IF(VALUE.GT.HIGH) THEN 

WRITE (TOUT, 30) NAME, HIGH 
CALL DELAY(6) 

GO TO 100 
ELSE 
END IF 
END IF 
RETURN 
END 
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CDAT85 



•k k 

* SUBROUTINE CDAT85 (LOW, HIGH .VALUE , NAME , IMOD) * 

* REAL VERSION * 

* HANDLE INPUT OF PARAMETER DATA FOR THE PROGRAM. INCLUDES SOME * 

* ERROR CHECKING. * 

* * 

* VARIABLES USED ARE: * 

* LOW : MINIMUM VALUE THAT WILL BE ACCEPTED, INPUT * 

* HIGH : MAXIMUM VALUE THAT WILL BE ACCEPTED, INPUT * 

* VALUE : VALUE READ FROM SCREEN, OUTPUT * 

* NAME : SCREEN PROMPT STRING, INPUT * 

* IMOD : FLAG USED IN ANOTHER PROGRAM, SET TO 1 ANYTIME * 

* CDAT85 IS CALLED., OUTPUT * 

* ESC : ESCAPE CHARACTER, (27) * 

* * 
•k-k-k-k'k-k-k^k-k'k-k-k-k'k-k-k-k^^ic-k-k'k-k-kic-k-k-kic-k'k-k-k-k'k-k-kic-kic'k-fc-k-k-k-k-k-k-k'k-k-k-k-k'k-k-kic'k-k-k-k'k-k'J^'k-k'kis-k 

SUBROUTINE CDAT85 (LOW , HIGH , VALUE , NAME , IMOD) 

COMMON/FI LS/IIN , lOUT , JOUT , KOUT , LOUT , MOUT , NOUT 
REAL LOW 

CHARACTER* 3 5 NAME 
CHARACTER*! ESC 
ESC = CHAR(27) 

10 FORMAT(' ENTER THE ' ,A35) 

20 FORMAT(' THE MINIMUM VALUE FOR ',A35,/' IS ',F6.2) 
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30 FORMATC THE MAXIMUM VALUE FOR ',A35,/' IS '.F6.2) 



* 








* 


•k 


CLEAR SCREEN, 


POSITION CURSOR, 


AND 


* 


k 


SET IMOD FLAG 


AND WRITE PROMPT 


TO SCREEN. 




k 








* 



*************************************************>V**************'*'*'***'*'* 
100 WRITE (lOUT,*) ESC,'[2J' 

WRITE (lOUT,*) ESC,'[10;0H' 

IMOD = 1 

WRITE (lOUT, 10) NAME 

****************************************■*•***■*•***■*•■*■**■*•***■********•!!:*****■*• 



* * 

* READ SCREEN VALUE. ON ERROR OR END OF DATA SET, GO BACK TO * 

* PROMPT . * 

* * 



*********************************************************************** 
READ(IIN,*,ERR = 100, END = 100) VALUE 



•jSr * 

* IF VALUES ARE OUT OF RANGE , PRINT ERROR MSG AND RETURN TO * 

* PROMPT . * 

* * 



IF (VALUE.lt. LOW) THEN 

WRITE (IOUT.20) NAME, LOW 
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CALL DELAY (6) 

GO TO 100 
ELSE 

IF(VALUE.GT.HIGH) THEN 

WRITE (IOUT.30) NAME, HIGH 
CALL DELAY (6) 

GO TO 100 
ELSE 
END IF 
END IF 
RETURN 
END 
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GAUSS 






ic ■sSr 

* SUBROUTINE GAUSS (NRHS) * 

* * 

* SOLUTION OF LINEAR ALGEGRAIC SYSTEM BY GAUSS ELIMINATION * 

* WITH PARTIAL PIVOTING. TAKEN FROM 'AN INTRODUCTION TO * 

* THEORETICAL AND COMPUTATIONAL AERODYNAMICS', JACK MORAN, * 

* JOHN WILEY & SONS, NEW YORK, 1984, PG 78. * 

* * 

* VARIABLES USED ARE: * 

* A() : COEFFICIENT MATRIX, AUGMENTED WITH RIGHT HAND * 

* SIDES IN COLUMN NEQNS+1, INPUT/OUTPUT , UNUSABLE * 

* AFTER RETURN * 

* NEQNS : NUMBER OF EQUATIONS IN MATRIX, INPUT * 

* NRHS : NUMBER OF RIGHT HAND SIDES, INPUT * 

* * 



******************•***•*■•***•*•*■*•*•*■■*■**■********■***********•***■*■*******•*******■* 
SUBROUTINE GAUSS (NRHS) 

COMMON/COF/A (10 1,101) , NEQNS 
NP = NEQNS+1 
NTOT - NEQNS+NRHS 

* * 

* SEARCH FOR THE LARGEST ENTRY IN THE (I-l)TH COLUMN, ON OR * 

* BELOW THE MAIN DIAGONAL. * 
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•:k"k:k':ki^ic'k-k'kic':k-k-k-:k-:k‘k‘:ki('k':k'k'k:k'k-k'ki(:kic'k:k’k-k'kic:k’k:k'k’k:k-k:k'k’k'k’k:kic:kic:ki(:'k':k'k'k'kic'kic^ 



DO 150 I = 2.NEQNS 
IM = I-l 
IMAX = IM 

AMAX = ABS(A(IM,IM)) 

DO 110 J = I.NEQNS 

IF (AMAX.GE.ABS(A(J,IM))) GO TO 110 
IMAX = J 

AMAX = ABS(A(J,IM)) 

110 CONTINUE 

'k'k'kic'k'k'k'k'k-k'kicic'ki('k'ki(’ki:'k':k:k'k'9c’:kic-:k:k'k’Jc’k’:k':k'k':k'k:k'k'k'k-:k'Jc':kic':k’k':k'kic'k^ 

* * 
* SWITCH THE (I-l)TH AND IMAXTH EQUATIONS * 



*********************************************************************** 
IF (IMAX.NE.IM) GO TO 140 
DO 130 J = IM.NTOT 
TEMP = A(IM,J) 

A(IM,J) = A(IMAX,J) 

A(IMAX,J) = TEMP 
130 CONTINUE 

-k':k'kici(ic'k-:k‘k'k-kic-k-k'k'k'k-k-k-kic'kic'k'k-kic-^c'k-k'kic-k'k'^c-k-k-k-k-k-k-k-kic-k-k'k'k'kic'k'k'k'k'^ 

ic * 

* ELIMINATE THE (I-l)TH UNKNOWN FROM ITH THROUGH (NEQNS)TH * 

* EQUATIONS . * 



140 DO 150 J = I.NEQNS 

R = A(J,IM)/A(IM,IM) 

DO 150 K = I.NTOT 
A(J,K) = A(J,K)-R*A(IM,K) 

150 CONTINUE 

^ic'kic'kic’:k-:k-:k-k'ki('k:k-k-:kic'k-:kic'k^-:ki('k-:k'k-:k-kic-k'kic-^ic-:k-k-k:k-:k-ki('k‘k'kic-k-:kicic-k'k'k'k'k'k'k'kic'k'k'kic'k 
'k 'k 

* BACK SUBSTITUTE * 



icic'k'kic'k'kic'kiicic'kicicicic'k'k'kicic'kic'k'kic'k'k'kic'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

DO 220 K = NP.NTOT 

A(NEQNS,K) “ A(NEQNS,K)/A(NEQNS,NEQNS) 

DO 210 L = 2.NEQNS 
I = NEQNS+l-L 
IP = I+l 

DO 200 J = IP.NEQNS 

A(I,K) = A(I,K) - A(I,J)*A(J,K) 

200 CONTINUE 

A(I,K) = A(I,K)/A(I,I) 

210 CONTINUE 

220 CONTINUE 
RETURN 
END 
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RNDUP 



'kii:i^ii:'k-k'k:k'k'k'k'k:k'k'k'k:k'k'ki^'ki^'kic'k:k:k'k:k'k:k^'k'k'k'k'k'k'k'k:kic'k'k'k'k^'k-k:k-k'k'k'k’k'k'k'k-k'k'k-k'k-k-k'k'k 



"k k 

* SUBROUTINE TO ROUND UP THE VALUE OF AN AXIS LIMIT. * 

^ FOR EXAMPLE IF THE VALUE PASSED IS .0018, THE VALUE RETURNED ^ 

^ IS ,0020. THIS IS NECESSARY TO ENSURE THAT AXIS LIMITS ARE * 

^ WELL DEFINED. 

k k 



k 

k 

k 

k 

k 

k 

k 

k 

k 

k 



VARIABLES USED ARE: ^ 

VALIN THE VALUE PASSED ^ 

VALOUT THE VALUE RETURNED ^ 

VAL A TEMPORARY VARIABLE * 

IRTN AN INTEGER MULTIPLIER USED TO RETURN TO THE CORRECT ^ 

ORDER OF MAGNITUDE. ^ 

IS IDE AN INTEGER ADDED TO ROUND UP. DEPENDS ON WHETHER THE ^ 

VALUE PASSED IS POSITIVE OR NEGATIVE. * 

IV A TEMPORARY VARIABLE ^ 

k 



kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

SUBROUTINE RNDUP (VALIN, VALOUT) 

RTN =1.0 
VAL = VALIN 
IF (VAL. EQ. 0.0) THEN 
VALOUT = VAL 
GOTO 200 
ELSE 
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END IF 



IF (VAL.LT.0.0) THEN 

ISIDE = -1 
ELSE 

ISIDE = 1 
ENDIF 

100 IF (ABS(VAL) .GE.10.0) THEN 

VAL = VAL/10.0 

RTN = RTN*10.0 
ELSE 

IF(ABS(VAL) .LT.1.0) THEN 
VAL = VAL*10.0 
RTN = RTN/10.0 

ELSE 

IV = AINT(VAL)+ISIDE 
VALOUT = REAL(IV*RTN) 
GO TO 200 

ENDIF 
ENDIF 
GO TO 100 
200 RETURN 
END 
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SDIG 



'k'k-krk'k-kic-k:k-k:kic'k:k-k:k:kic’:k’k’k’:kic-:k'k'k:k'k:k'k:k:k'k'k'k'k:k'k:k'k'k'ki^'k'k'k:k'k:k:k:k:k'k:k'kicic'k’:fcic'k:kic 



•A- 






* 


•A- 


SUBROUTINE SDIG (MININ , MINOUT , MAXIN , MAXOUT , IDIG) 


k 


ic 






k 


* 


SUBROUTINE TO GET THE RANGE OF VALUES IN THE AXIS LIMITS 


k 


•k 


TO BE 5 


UNITS APART. AND RETURN THE NUMBER OF SIGNIFICANT 


k 


k 


DIGITS 


FOR THE AXIS. 


k 


k 






k 


k 


VARIABLES USED ARE: 


k 


k 


MININ 


THE MINIMUM AXIS VALUE PASSED 


k 


k 


MINOUT 


THE MINIMUM AXIS VALUE RETURNED 


k 


k 


MAXIN 


THE MAXIMUM AXIS VALUE PASSED 


k 


k 


MAXOUT 


THE MAXIMUM AXIS VALUE RETURNED 


k 


k 


VAL 


A TEMPORARY VARIABLE 


k 


k 


VALNOT 


A TEMPORARY VARIABLE 


k 


k 


RTN 


A MULTIPLIER USED TO RETURN TO THE CORRECT 


k 


k 




ORDER OF MAGNITUDE. 


k 


k 


ISIDE 


AN INTEGER TO DETERMINE WHETHER THE VALUE IS POSITIVE 


k 


k 




OR NEGATIVE 


k 


k 


IV 


A TEMPORARY VARIABLE 


k 


k 






k 



':k'kicic'kicrk-kic'k'kic'kic'kic-k-k‘k'k-kic-kic-kicicic'k'k'k-kicic'kii:'kicic':kic-k'k'kic'k'k'k'kicicicic-k 

SUBROUTINE SDIG (MININ , MINOUT , MAXIN , MAXOUT . IDIG) 
REAL MININ , MINOUT , MAXIN , MAXOUT , RTN , VAL , VALNOT 
INTEGER IDIG , ICASE, ISIDE , IV , ITEST 
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RTN =1.0 

IF( ABS (MININ) .LT.ABS(MAXIN)) THEN 

* * 

* ICASE = 1 IS FOR THE MIN AXIS VALUE BEING SMALLER IN MAGNITUDE * 

* THAN THE MAX AXIS VALUE. * 

* * 
A*****************************************************'*''*'**'*'**'***'*'*'**'*'** 

ICASE = 1 
VAL = MAXIN 
VALNOT = MININ 
ELSE 

•k-k-k-k-k-k-k-k-k-k-kidrk'k'kickirkic^-k-k-kic-k'k-k-k-kidcic'kick-kirk-k-k'kidrJcidck'k-k-k-k'k'k'k-k-k'k-kirk-k-k'k'k-k-kic'k'k 
* * 

* ICASE = 2 IS FOR THE MAX AXIS VALUE BEING SMALLER IN MAGNITUDE * 

* THAN THE MIN AXIS VALUE. * 

* -k 
•kk'kk-kkkk-kkkkkkkkk-kkkk-kkkk-kkkkkkkkkkkkkkk-kkkkkkkkkkkkkkkkkkkkkk-k-k-kk-kkk'k-k 

ICASE = 2 
VAL = MININ 
VALNOT = MAXIN 
ENDIF 

kkkkk-kkirkk-kirk-kk-kkkk-kkk-k-k-kk-k-kkkkkkk-k-kkkkkirkkk-k-k-k-k-k-k-k-kk-k-k-kk-k-k-k-k-kkk-k-k-kkkk-k 



DETERMINE WHETHER THE SMALLER (IN MAGNITUDE) AXIS VALUE IS 
POSITIVE OR NEGATIVE. 
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ISIDE = SIGN(l.VALNOT) 

100 IF (ABS(VAL) .GE.10.0) THEN 
VAL = VAL/10.0 
RTN = RTN*10.0 
ELSE 

IF(ABS(VAL) .LT.1.0) THEN 
VAL = VAL*10.0 
RTN - RTN/10.0 

ELSE 

IF (VALNOT.EQ.0.0) GOTO 300 
IV = NINT(VALNOT/RTN) 

IF(IV.EQ.O) IV - IV + ISIDE 
GO TO 200 

END IF 
END IF 
GO TO 100 
200 CONTINUE 

ITEST = ABS(NINT(VAL) -IV) 

IF ((ITEST. EQ. 5) .OR. (ITEST . EQ . 10) .OR. (ITEST. EQ. 15) .OR. 

* (ITEST.EQ.20) .OR. (NINT(VAL) .EQ. -IV)) GOTO 300 
IV - IV + ISIDE 

ITEST = ABS(NINT(VAL)-IV) 

IF ((ITEST. EQ. 5) .OR. (ITEST. EQ. 10) .OR. (ITEST . EQ . 15 ) .OR. 

* (ITEST.EQ.20) .OR. (NINT (VAL) .EQ. -IV)) GOTO 300 



VAL = VAL - ISIDE 



ITEST = ABS(NINT(VAL)-IV) 

IF ((ITEST.EQ.5) .OR. (ITEST . EQ. 10) .OR. (ITEST . EQ . 15) .OR 

* (ITEST.EQ.20) .OR. (NINT(VAL) .EQ. -IV)) GOTO 300 
IV = IV -(- ISIDE 

ITEST = ABS(NINT(VAL) -IV) 

IF ((ITEST.EQ.5) .OR. (ITEST. EQ. 10) .OR. (ITEST. EQ . 15) . OR 

* (ITEST.EQ.20) .OR. (NINT( VAL) .EQ. -IV)) GOTO 300 
VAL = VAL - ISIDE 

300 IDIG = NINT(LOG10(1.0/RTN)) + 1 
IF (ICASE.EQ.l) THEN 
MINOUT = IV*RTN 
MAXOUT = NINT(VAL)*RTN 
ELSE 

MAXOUT = IV*RTN 
MINOUT = NINT(VAL)*RTN 
ENDIF 
RETURN 
END 
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APPENDIX B. NOMENCLATURE 



AR aspect ratio 

b span 

c local chord 

c average chord 

c ^ root chord 

c tip chord 

Cj) citag coefficient 

induced drag coefficient 
C]^ lift coefficient 

C]^l ideal lift coefficient 

moment coefficient about the leading edge of an 
airfoil 

Cm root chord for circulation model 

moment coefficient 

^Mo moment coefficient about the y axis 

^Mac moment coefficient about the aerodynamic center 

AC p pressure difference coefficient 

d gtid element inset 

D dt ag 

AD I section induced drag 

AF force per unit span 

K vortex drag factor 

L lift 
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M 


moment 


r 


distance between field and control points 


S 


planform area 


^eff 


local effective velocity 


^co 


remote velocity 


w 


induced velocity 


X 


chordwise coordinate 


^ac 


X coordinate of aerodynamic center 


^cp 


chordwise position of center of pressure 


y 


spanwise coordinate 


z 


height coordinate 


(x,y) 


field point 


(Xo^Yo) 


control po int 


a 


angle of attack 


r 


circulation 


AT 


incremental circulation 


A 


leading edge sweep angle 


e 


transformed spanwise coordinate 


X 


tangent of leading edge sweep angle in circulation 
model 


X 


taper ratio in VORTEX model 


p 


dens i ty 


T 


percent wing taper in circulation model 


4> 


transformed chordwise coordinate 
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control point 
field point 
wing shape 



position where velocity is evaluated 
position of the vortex element 
surface slope of the wing, including twist 
from root to tip and camber 
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